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ABSTRACT 


Optimization algorithms were developed for use with the Adaptive Joint Time- 
Frequency (AJFT) algorithm to reduce Inverse Synthetic Aperture Radar (SAR) image 
blurring caused by higher-order target motion. A specific optimization was then applied 
to 3D motion detection. Evolutionary search methods based on the Genetic Algorithm 
(GA) and the Particle Swarm Optimization (PSO) algorithm were designed to rapidly 
traverse the solution space in order to find the parameters that would bring the ISAR 1m- 
age into focus in the cross-range. 3D motion detection was achieved by using the AJTF 
PSO to extract the phases of 3 different point scatterers in the target data and measuring 
their linearity when compared to an ideal phase for the imaging interval under investiga- 


tion. The algorithms were tested against both simulated and real ISAR data sets. 
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EXECUTIVE SUMMARY 


Inverse synthetic aperture radar ISAR) imaging 1s an effective way to acquire 
high resolution images of targets of interest at long range and as such is an irreplaceable 
tool in the task of non-cooperative target recognition (NCTR) of both ships and aircraft. 
The Canadian Air Force is currently upgrading her fleet of CP-140 Aurora maritime pa- 
trol aircraft to possess NCTR through ISAR 1n order to increase the capability of the Ca- 
nadian Forces in both sovereignty patrols of Canadian territory and the protection of the 
Canadian Patrol Frigates and allied ships operating abroad as a coalition. The United 
States Navy’s equivalent aircraft, the P-3 Orion, currently possesses this capability. 
Therefore, effective ISAR imaging will have a real impact in the decision-making proc- 


ess in future military operations involving both Canadian and American forces. 


One significant problem with ISAR image formation is the assumption of time in- 
variance of the Doppler frequency used to resolve the image in the cross-range. Time- 
variant Doppler becomes present in an ISAR signal when an aircraft is maneuvering or a 
ship is pitching and rolling during the coherent processing interval and is typically re- 
ferred to as motion error. Because of this motion error and the fast Fourier transform’s 
(FFT) basic assumption of stationary signals, the use of the FFT to resolve the image in 
the cross-range will cause extensive blurring and leave the image unrecognizable even to 


the most experienced ISAR operators. 


One proven method of motion compensation is the adaptive joint time-frequency 
(AJTF) algorithm. However, this algorithm is not without two significant weaknesses of 
its own. The first problem is that the computational burden of the exhaustive search used 
to extract the motion compensation parameters is quite large which limits its usefulness 
in an operational situation. The second problem with the AJTF algorithm 1s that the tar- 
get must follow the mathematical model of the AJTF which, in particular, states that rota- 
tional motion must be confined to a 2D plane (1.e., no 3D motion). 

To address the first problem, two evolutionary searches, a Genetic Algorithm 


(GA) and a Particle Swarm Optimization (PSO), were investigated, developed and evalu- 


X1X 


ated in order to take advantage of their properties of rapid convergence. Both the GA and 
the PSO performed excellently during their testing against simulated and real radar data 
in the extraction of the search parameters from the solution space. They would typically 
converge to the search parameters required for focusing in less than 50% of the cost (or 
fitness) function calculations required by the exhaustive search. Also, the PSO proved to 
be a much more efficient algorithm for this optimization problem with its performance 


consistently exceeding that of the GA. 


To address the second problem, the PSO algorithm designed above was combined 
with the 3D motion detection algorithm developed at DRDC — Ottawa by Thayaparan. 
The resulting algorithm was able to effectively and efficiently sort good imaging intervals 
in the real radar data set from poor imaging intervals which violated the 2D motion as- 
sumption of the mathematic model. It was shown that the imaging interval indicated as 
possessing a low degree of 3D motion created well-focused images after being motion 
compensated. Imaging intervals that were indicated as possessing a high degree of 3D 
motion did not focus even after being motion compensated. These imaging intervals can 


now be discarded instead of presented to the ISAR operator. 


This thesis draws from many sources 1n the radar signal processing community 
and cites them appropriately. There are two contributions made by this thesis that are 
novel in their entirety. First, the automatic detection of focus using the FFT for determin- 
ing basis function suitability in the AJTF algorithm is new, effective and fast and not 
used 1n any of the cited references. Secondly, the AJTF PSO algorithm, with the PSO 
based on the work in [7], 1s unique to its application to ISAR image focusing. The work 
on 3D motion detection draws heavily from the work 1n [12] and provides two significant 
improvements. First, the measure of non-linearity in an imaging interval is based on the 
percentage of deviation from the line of least-squares fit as opposed to the basing it on the 
difference between the line of least-squares fit and the actual line depicting the linearity 
of phases. Secondly, the AJTF PSO developed in the first part of this thesis 1s applied to 
the 3D motion detection algorithm. This creates an algorithm that evaluates the degree of 
3D motion in an imaging interval both rapidly and accurately which 1s useable by an 
ISAR operator in an operational situation. 
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I. INTRODUCTION 


A. RESEARCH REQUIREMENT AND OBJECTIVES 

Inverse synthetic aperture radar ISAR) imaging 1s an effective way to acquire 
high resolution images of targets of interest at long range and as such 1s an irreplaceable 
tool in the task of non-cooperative target recognition (NCTR) of both ships and aircraft. 
The Canadian Air Force is currently upgrading her fleet of CP-140 Aurora maritime pa- 
trol aircraft to possess NCTR through ISAR 1n order to increase the capability of the Ca- 
nadian Forces in both sovereignty patrols of Canadian territory and the protection of the 
Canadian Patrol Frigates and allied ships operating abroad as a coalition. The United 
States Navy’s equivalent aircraft, the P-3 Orion, currently possesses this capability. 
Therefore, effective ISAR imaging will have a real impact in the decision making process 


in future military operations involving both Canadian and American forces. 


As explained in [1], [2] and [3] and repeated in Chapter II of this thesis, one 
drawback of ISAR is that higher order target motion, observed when ships are pitching 
and rolling in rough seas or when aircraft are maneuvering, introduces time-varying Dop- 
pler into the radar signal. It is this time-varying Doppler that causes severe image blur- 
ring in the cross-range during the ISAR signal processing stage and effectively renders 
the target’s image unrecognizable even to the most experienced ISAR operators. There- 
fore, for ISAR to realize the operational goal of NCTR, methods must be developed to 
compensate for the higher order motion. The first objective of this thesis was to continue 
the work by Thayaparan [1] and Chen [3] on the Adaptive Joint Time — Frequency Algo- 
rithm (AJTF) and to examine different methods of efficiently removing time-varying 
Doppler from the radar signal in order to form images that are well focused in the cross- 


range. 


The AJTF algorithm in [1] and [3] assumes 2D motion on a rotational axis during 
the required ISAR dwell interval. Unfortunately, true targets do not always follow this 
assumption and within the radar signal there exists 3D motion which cannot be compen- 
sated by the AJTF algorithm. One way to deal with the existence of 3D motion in an 


ISAR imaging data set is to determine to what degree 3D motion exists in an imaging 1n- 


terval of the data set. The imaging intervals with a small degree of 3D motion are imaged 
while the intervals that possess a high degree of 3D motion are discarded. The research 1n 
[3], [4] and [12] provide the background for 3D motion detection using the nonlinearity 
of phase method. Therefore, the second objective of this thesis is to apply the fast AJTF 
algorithm found in the first part of this thesis to 3D motion detection methods found in 
[3], [4] and [12]. 
B. SCOPE AND ORGANIZATION 

1. Scope 

This thesis uses two ISAR data sets, one which 1s a B727 simulated data set cre- 
ated by Chen [18] and 1s designed to display the ideal motion error found 1n a blurred 
ISAR image and is available from the Naval Research Laboratory website. An example 
of the unfocused and focused image is found in Figure 1. The other ISAR data set, which 
was used 1n [1], is the 6 — Point Delta Wing experimental data set collected by Defence 
Research and Development Canada — Ottawa (DRDC — Ottawa) at their radar range lo- 
cated at Shirley’s Bay, Ontario. A photo of the 6 — Point Delta Wing is shown in Figure 


2 and an example of its blurred and focused image can be found in Figure 3. 
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Figure 1. B727 Simulated Data Set - Unfocused and Focused Images. 
Z 





Figure 2. _ Picture of 6 — Point Delta Wing (From [1].) 
Note that one reflector is offset to give a sense of orientation. 
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Figure 3. | Example of an Unfocused and Focused 6 — Point Delta Wing. 
These two data sets were used to test the algorithms designed in this thesis. The 


B727 simulated data set is an excellent data set since its ideal motion error characteristics 
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make it perfect for testing the first cut of an imaging algorithm. The 6 — Point Delta 
Wing data set is an experimental data set and it 1s particularly useful since the entire data 
set consists of 60,000 pulses in the cross-range. The 60,000 pulses can be cut into differ- 
ent size imaging intervals with each of these intervals displaying a different amount of er- 
ror. Also the 6 — Point Delta Wing experimental data set is the only one of the two data 
sets that can be used for the 3D motion detection portion of the thesis. 

Zz: Primary Research Questions 

Based on the objectives above, there are two research questions that will be the 


primary subject of this thesis: 


a) Can the effective application of a genetic algorithm or a particle swarm 
optimization algorithm on the parameter search of the basis function for the AJTF algo- 


rithm significantly reduce the computation burden required to focus the ISAR image? 


b) Can the fast AJTF algorithm found in the first part of this thesis be applied 
to the 3D motion detection problem in order to rapidly sort the good imaging intervals, 
which obey the 2D motion assumption, from the imaging intervals that contain a high de- 
gree of 3D motion and cannot be focused? 

3: Organization 


This thesis is organized into six parts. 


iB Chapter II introduces ISAR signal processing to the extent required to un- 


derstand the Adaptive Joint Time-Frequency Algorithm (AJTF); 


il. Chapter HI discusses time-frequency transforms and the AJTF algorithm, 


as found in [1], [2] and [3]; 


111. Chapter IV presents the Genetic Algorithm (GA) and the Particle Swarm 
Optimization (PSO) Algorithm that were written to optimize the AJTF al- 


gorithm along with the ISAR imaging results; 


IV. Chapter V presents the algorithm for 3D motion detection and ISAR imag- 
ing results of imaging intervals determined to possess a high and low de- 


gree of 3D motion. 


V. Chapter VI summarizes the results that were achieved 1n this thesis and of- 


fers conclusions and recommendations for future work; and 


V1. The Appendix lists the MatLab code that was written for this thesis. 
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H. INTRODUCTION TO INVERSE SYNTHETIC APERTURE 
RADAR (ISAR) PROCESSING AND THE FOCUSING PROBLEM 


This chapter introduces the basic equations required to understand ISAR imaging 
and the focusing problem, which occurs when non-stationary signals such as time- 
varying Doppler are present in the ISAR return signal. It shows how Doppler resolution 
is related to the coherent processing interval (CPI) and how this relates to the amount of 


motion error in an image. 


For the purposes of this thesis, the ISAR return signal 1s at the point in the proc- 
essing where all of the point scatterers have been placed in their proper range bins and a 
fast Fourier transform (FFT) will focus the scatterer in their appropriate cross-range bins. 


The ISAR signal without error can be represented as [1,3]: 





2E Jo | Ry +X, cos(@(t)) + y, sin(0(0))]}. (1) 


Nx 
s(x,t)= A exp |=) 

k=l C 
where x is the range bin, NV, is the number of scatterers in the range bin x, A, is the mag- 
nitude of the radar return signal, R, 1s the distance to the rotational center of the target 
which is constant during the coherent processing interval, Ee y,) is the down-range and 
cross-range distance from the target’s rotational center to the k” point scatterer. Since no 
motion error is being assumed at this point, @(¢) can be written as [1], [3]: 


O(t) =9, + M1, (2) 


where ©2 is the constant angular velocity of the radar target from which the Doppler 


cross-range location of the scatterer is extracted. If O(¢) is assumed to be small, Equa- 


tion (1) reduces to [3]: 
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which can easily be focused by the fast Fourier transform since the extracted Doppler 1s 
time-invariant. The time-invariance of the Doppler is critical for successful focusing us- 


ing the fast Fourier transform. 


The second equation of importance deals with the Doppler resolution, Af,,, and 
its associated effect on the cross-range resolution, Ar... Doppler resolution is associated 
with the coherent processing interval, 7,,,, such that [3]: 

Af, Dan (4) 


7, 


CPI 


which states that as the coherent process interval is increased, the Doppler resolution in- 


creases. Similarly, resolution in cross-range can be defined as [3]: 


Re a 5 
(ate = sag a ” 


Therefore, from Equation (5), it is clear that, as the coherent processing interval 1s 1n- 








creased, the size of the Doppler bins 1n the cross-range becomes smaller and it is possible 
to resolve two different point scatterers in the same range bin even if they are only 
slightly separated in Doppler. In this thesis, the number of pulses in the cross-range 1s 


analogous to 7,,,. A greater T..,, leads to a greater number of pulses in the cross-range 


and, subsequently, an ISAR signal with greater Doppler separation between point scatter- 


ers that are collocated in the same range bin. 


If Equation (3) held for real ISAR returns from fast moving targets, 1t would make 


sense to use the maximum number of radar pulses (by maximizing the 7,,, ) to form the 
high resolution image. However if the entire 60,000 pulses ( 7..,, = 30 seconds) of the 6 — 


Point Delta Wing experimental data set is focused using the fast Fourier transform, the 
resulting image, as shown in Figure 4, is extremely blurred and bears no resemblance to 


the actual target in Figure 2. Even if a more reasonable 7,,, 


is chosen, such as a T7.,, = | 
second with 2048 pulses in the cross-range, the motion error is still significant enough to 
cause the image to blur beyond recognition. This can be seen by referring to the unfo- 


cused image in Figure 3. 


If, instead, the 7 


cp, 18 reduced to 0.128 seconds so that there are now only 256 
pulses in the cross-range, motion error can be avoided. The side effect of this choice 1s 


that the cross-range resolution is greatly reduced. Most imaging intervals will generate 
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images that resemble those found in Figure 5. These images possess no noticeable error 
in the cross-range but also possess no resolution in the cross-range with most point scat- 


terers collapsed together on top of each other in the center Doppler bin. 


6- Point Deta Wing Data Set 
£0 000 Pilices in (noe. Range - Tip = 30 seconds 


Down Range - Range Bins 





2.94 2.96 2.98 3 3.02 3.04 3.06 
Cross Range - Doppler Bins + 10" 


Figure 4. 60,000 Pulse 6 — Point Delta Wing Data Set Focused with the Fast Fou- 
rier Transform. 
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Figure 5. 256 Pulse Imaging Interval 6 — Point Delta Wing Data Set Focused with 
the Fast Fourier Transform. 
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Therefore, our motivation for a fast and effective motion compensation algorithm 


becomes obvious. Any imaging interval with a 7/,, that leads to good separation in Dop- 


pler of the point scatterers within the same range bin will possess excessive motion error, 


time-varying Doppler and a blurred image. An imaging interval with a small 7,,, will 


have no error, time-invariant Doppler but also no separation of point scatterers in the 
cross-range. What is required is an algorithm that can remove the time-varying Doppler 


from an imaging interval of sufficiently long 7.,, such that the resulting image is both 


well focused and well separated in Doppler. 


With our requirement for an algorithm to correct motion error in an ISAR signal 
now clear, we can proceed to Chapter III where time-frequency transforms and the Adap- 
tive Joint Time-Frequency algorithm are introduced. Time-frequency transforms and the 
AJTF are both crucial tools in achieving effective motion compensation in an ISAR sig- 


nal. 
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Hl. INTRODUCTION TO TIME-FREQUENCY TRANSFORMS 
AND THE ADAPTIVE JOINT TIME-FREQUENCY ALGORITHM 


As stated in Chapter II, in order to form good ISAR images that are well resolved 


in the cross-range, there is a requirement for a coherent processing interval, 7_,,, that is 


Pl? 


sufficiently long such that the cross-range resolution, Ar,., is acceptably small. Unfortu- 


nately a T.,, that leads to a small Arv.. also leads to blurring in the cross-range since the 


Doppler of the target’s scatterers 1s no longer stationary, which 1s a critical requirement 
of the FFT. In this chapter, time-frequency transforms are introduced so that the non- 
stationary nature of the time-variant Doppler can be understood. Also, the AJTF algo- 
rithm, which 1s an effective way to compensate for time-variant Doppler is introduced. 
A. TIME-FREQUENCY TRANSFORMS 

As stated in [1], [3] and [5], the one weakness of the fast Fourier transform (FFT) 
is that it loses all information on how the signal changes with time. As explained in 
Chapter I, ISAR signal processing that uses the assumption of time-invariant Doppler 


during the coherent processing interval, 7..,,, will inevitably introduce severe blurring 


CPI? 
into the final high resolution image. Therefore, it becomes necessary to introduce the 
time-frequency transform and how it can be used to create an analogy between what is 
occurring in the time-frequency plane, the Fourier Spectrum and the radar image. 

1. Short Time Fourier Transform 

As discussed in [5], the short time Fourier transform (STFT) is a unique and 
straightforward way of introducing a time dependency into the Fourier transform. This 


is done by pre-windowing the FFT around a certain time instance and can be defined as 


[Si 


S.(t.»,h)= [ s.(Wh (ue Pou (6) 


—0O 


where S_ is the STFT of range bin x, s, is the ISAR signal for range bin x and h’ (y— 1) is 


the complex conjugate of the windowing function, which for purposes of this thesis is the 


default Hamming window. Figure 6 shows the time-domain waveform of an arbitrary 
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ISAR radar signal and the Hamming window. To generate the STFT, the Hamming win- 
dow is moved across the radar signal and at each time instance the FFT is taken of the re- 


sult of the multiplication of the window with the ISAR radar signal. 


Using the time-frequency toolbox from [5], the STFT of an ISAR radar signal can 
be generated and the analogy between what is shown by the STFT, seen in the power 
spectrum of an ISAR range bin and displayed in the radar image, can be illustrated. Fig- 
ure 7 shows these three graphs when taken from a blurred ISAR image. By examining 
range bin 31 in the close-up of the B727 blurred ISAR image (left), it can be seen that 
there appears to be two scatterers smeared across several Doppler bins in the cross-range. 
The same deduction can be made by examining the power spectrum (upper right), as it 
can be seen that there are two Doppler smeared point scatterers located in this range bin. 
Finally, the STFT (lower right) can be examined. From the STFT, it can be very clearly 
seen that there are two point scatterers located in range bin 31. The STFT also shows the 
nature of the time-varying Doppler of the two point scatterers as the ISAR signal is inte- 
grated from pulse to pulse during the coherent processing interval. The first point scat- 
terer’s Doppler decreases from pulse to pulse while the second point scatterer’s Doppler 


increases from pulse to pulse. 


Figure 8 shows the same information as Figure 7; however, it is for a focused 
ISAR image. Examining the close-up ISAR image from a focused B727 data set, it can 
be seen that the point scatterers in range bin 31 have collapsed in the cross-range and are 
well focused. It can also be seen that these scatterers have been shifted in Doppler. 
Since this shift is applied equally across the image, it has no adverse affect. More impor- 
tant are the results that are seen in both the STFT and the power spectrum of range bin 
31. The STFT shows that the two point scatterers, which are now lines with zero slope in 
the time-frequency plane, no longer possessing time-varying Doppler during the coherent 
processing interval. This characteristic is echoed in the power spectra of the two point 
scatterers which now come to sharp points rather than exhibiting smearing of their power 


spectra across several Doppler bins as occurred in the unfocused image. 
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Figure 6. _ISAR signal, s,., of Range Bin 31 and the Hamming Window (After 
[5].) 
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Figure 7. Blurred ISAR Image (close-up of Range Bin 31 in B727 data set), Power 
Spectrum of Range Bin 31 and its STFT. 
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Figure 8. | Focused ISAR Image (close-up of Range Bin 31 in B727 data set), 
Power Spectrum of Range Bin 31 and its STFT. 


In concluding this section, it 1s important to reiterate the following points as they 


will become critical in chapter 4 when image focusing is discussed: 


a) The range bins of a blurred ISAR image will have point scatterers that ap- 
pear as lines with slopes (or even curved for higher order motion) when transformed to 
the time-frequency plane and the slope 1s an indication of the presence of time-varying 
Doppler. This same range bin will have the peaks of its power spectrum, which are pro- 
duced by each point scatterer in the range bin, smeared across several Doppler bins when 


time-varying Doppler is present in the signal. 


b) The range bins of a focused ISAR image will have point scatterers that 
appear as lines with a slope of zero in the time-frequency plane. The power spectrum of 
these same range bins will have sharp, well defined peaks located solely in their appro- 


priate Doppler bins with minimum spillage into adjacent bins. 
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2. Weaknesses of the STFT and Time — Frequency Transform 

Despite the usefulness of the STFT and other time — frequency transforms, as 
listed in [5] and [3], there are two primary weaknesses which will be of concern in the fo- 
cusing problem of this thesis. The first weakness is the computational time of the STFT, 
which 1s the fastest time-frequency transform available in the time-frequency toolbox 
from [5]. Calculation time for the STFT is significantly longer than the FFT and this be- 
comes critical if the STFT must be computed multiple times. The second weakness of 
the STFT, as stated in [2] and [5], 1s that it possesses poor resolution and cannot be si- 
multaneously well resolved in both time and frequency. This is not a significant concern 
in the B727 simulated data set as all point scatterers in each range bin are well separated 
in Doppler. However, for the 6 — Point Delta Wing experimental data set, the point scat- 
terers are not well separated in Doppler and it 1s difficult to see how the Doppler of each 
point scatterer changes with time. To increase time-frequency resolution, as recom- 
mended 1n [3] and [5], a different time-frequency transform could be chosen; however, 
this typically results in an increase in computational time, or in the case of the Wigner- 
Ville Distribution (WVD), cross-terms appear in the time-frequency image. As shown in 
Figure 9, which is a WVD of range bin 31 of the unfocused B727 simulated data set, 
these cross-terms can make the interpretation of the time-frequency image very difficult 


and the automation of line detection close to impossible. 
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Figure 9. — Wigner-Ville Distribution of Unfocused Range Bin 31 (After [3] and 
[5]-) 
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B. ADAPTIVE JOINT TIME-FREQUENCY ALGORITHM 

The AJTF algorithm can be found in [1], [2] and [3] and is capable of removing 
time-varying Doppler in a radar signal that 1s due to translational and rotational motion 
error. As explained in the references, the capability of the AJTF is limited by the hold- 
ing of three primary assumptions. The first assumption 1s that the ISAR signal has been 
error corrected in the down-range direction and that each point scatterer is in its proper 
range bin. The second assumption is that of a rigid body, which implies that motion error 
that exists in one point scatterer will exist in all point scatterers of the target. This as- 
sumption 1s violated in the presence of jet engine modulation (JEM), helicopter rotors and 
the micro-Doppler phenomenon. The third assumption is that motion occurs in a 2D 
plane with respect to the radar. The AJTF algorithm is currently unable to correct for 3D 
motion. For this thesis, the first two assumptions will always hold. The assumption of 
2D motion does not hold for all imaging intervals of the 6 — Point Delta Wing experimen- 
tal data set and, as a result of the 3D motion error, these imaging intervals will not image 
well. As stated in the introduction, it is the second part of this thesis that deals with 3D 
motion error so, for now, it will be assumed that all three assumptions hold. 

1. The Motion Model for a Moving ISAR Target 

The motion model for an ISAR target is introduced in [2] and further explained in 
both [3] and [1]. The following explanation 1s drawn from all three sources but primarily 
[1] and [3]. Figure 10 shows the motion geometry of an ISAR target. Coordinates 


(U,V,W) is the radar coordinate system, while coordinates (x,y,z) is the coordinate 
system used to define the target and (X,Y, Z) is the coordinate system translated from 


the radar and is used to describe the target’s rotation. It 1s centered at the geometric cen- 


ter of the target, about which the target will rotate. If P(x, y) is the location of a point 


scatterer, then the initial distance from the radar to the k” point scatterer located on the 


target at t= 0 1s [3]: 


R, = R, +x, cos(@—a)+y,sin(@,-a), (7) 
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where R, is the initial line-of-sight distance from the radar to the target’s geometric cen- 
ter, G, is the initial angle between (Y,Y,Z) and the (x, y,z) coordinate system and @ is 


the azimuth angle of the scatterer off of the (U,V,W) axis. 
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Figure 10. Geometry of the Radar Imaging of a Target (From [3].) 


Now, if the target possesses a rotational motion, 02, about the Z axis and a trans- 


lational motion with radial velocity, V,, then both the range and rotational angle become 


a function of time and Equation (7) become time dependent. The time dependent 


R(t) can be expanded into a Taylor series polynomial as [1]: 


R(t)=R,+V pt + Vit + Vat tones (8) 





OV, ee ee ae ae 
where V, =——* =a, which is the acceleration of the scatterer in the radial direction and 
t 
the additional derivatives of V, are terms of higher order translational motion. The time 
dependent @(t) can be expanded into a Taylor series polynomial as [2]: 


7 
3! 


l 


2! 


O(t) =0,+Ot+—O't’ +—Q'"F +..., (9) 
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where 2)’ = = which is the initial angular acceleration of the point scatter and the addi- 


tional derivatives of Q are terms of higher order rotational motion. Taking Equations (8) 


and (9), substituting them into Equation (7) and arbitrarily setting @ and 6, to zero, we 


have [2]: 


l 


i: 1 
R,(t)=R, + Vat + Ve t+ ..4+X, cos[ 1+ 20 er 


Q"r _ 
; (10) 
+y, sin( r+ TOP +—O'P + ie 

2! 3! 


Now taking Equation (10) and calculating the Taylor series expansion of the sine and the 


cosine about ¢ yields [1]: 


1 
(Vi —Q?x,+Q%y, )P +... (11) 


2! 


R, (t) = (Ry + x, ) + (Qy, + Vz t+ 


which is the time-varying range to a one point scatterer on the target. 
2. The ISAR Signal and the Motion Compensation Algorithm 


Now Equation (11) can be substituted into Equation (1) which gives the equation 


for an ISAR return signal from a collection of NV, point scatterers located in range bin x 


where the ISAR return signal is not an ideal, errorless signal [1]: 


Nx 
saa) A exp |-/ AL (R, + x,)+(Qy, + V.)t+(V, -O?x, + Oy, )t? + Jf (12) 
k=] : 


If, as in [1], the following substitutions are made: 


Ay = Ky +X, 
2 =O: (13) 


] 
a, arr -O°x, +Q'y,), 


the following observations are made. First, a, defines the coordinate position of the scat- 


terer on the target. The second term, a,, 1s the radar line-of-sight translational motion 


and constant rotational motion and 1s linearly dependant on time. If all other terms are 


negligible, then a simple fast Fourier transform 1s sufficient to form a focused image. 
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The second term, a, (as well as subsequent terms if higher order motion is present), con- 
tains V, which is the translational motion error and 'y, which is the rotational motion 


error that has a dependency on the cross-range. The (Qr)’ can be neglected since [1] 
states that for centimeter and millimeter wave radars, this term is generally much less 


than 1. 


Therefore, we can define the phases of the uncompensated range bin as the phase 
function [1]: 
O(s)=a,+at+at +... (14) 
and since the phases 1n the range bin are time-varying, the instantaneous frequency can 


be expressed as [1]: 


1 0®(t) 
2x Ot 





f= (15) 


which is also time-varying. Therefore, the Doppler frequencies of the point scatterers in 


the range bin are time-varying and the image is blurred in the cross-range. 


The next step is to focus on a range bin and to find the translational motion com- 


pensation basis function of the form [3]: 
h, (t) = exp{—j2a O,(t)} 


| I. 2b. La (16) 
=exp)—j2a aegewe vat ae 


where ©,(¢) is the phase function for translational motion compensation and the parame- 
ters { fi1,f;>;3>---} are the AJTF search parameters for translational motion compensa- 


tion. A suitable basis function will closely resemble the prominent point scatterer in the 
time-frequency plane. Figure 11 shows an STFT of the unfocused range bin from the 

B727 simulated data set and the STFT of the best basis function for translational motion 
compensation (described later in this chapter). Once this basis function has been found, 
translational motion error can be corrected by applying the following equation from [1], 


[2] and [3]: 
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s(x,t), =s(x,t).*conjth, (1)} (17) 


where s(x,t), 1s the ISAR signal for all range bins with translational motion error re- 
moved, s(x,t) is the ISAR signal for all radar bins, .* defines array multiplication and 


conj 1s the MatLab complex conjugate operator. Figure 12 shows the same range bin in 
the time-frequency plane after translational motion compensation. Note that the top line, 
representing the prominent scatterer in the range bin has had its slope zeroed which 1ndi- 


cates the removal of time-varying Doppler. The remaining signal 1s now [1]: 





yee exp{ jtfotar, + 4y,0(0]| (18) 


C 


where O(t) 1s as defined in Equation (9) above. 





Figure 11. | Unfocused Range Bin (top) and the Best Basis Function /, (¢) for 
Translational Motion Compensation. 
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Figure 12. Time-Frequency Transform of Range Bin After Translational Motion 


Compensation. 


Now it is required that rotational motion compensation be completed on the sig- 


nal. This is done by choosing another range bin with a prominent point scatterer. If there 
is more than one scatterer in a range bin, it 1s possible to use the same range bin; how- 
ever, the prominent point scatterer used for translational motion must first be removed 
from the signal. The algorithms written for this thesis used separate range bins for trans- 
lational and rotational motion compensation. The first step 1s to again find a suitable ba- 


sis function as in Equation (16) and as shown in Figure 9: 


h, (t)=exp{—j2z-©,(t)}. (19) 


However, for rotational motion compensation, it is the phases of the basis function in 


Equation (19), ©,(¢), that are required [3]: 


l 1 
©,(t)= f,,t+ a fot’ + a fost? +... (20) 
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where © (tf) is the phase function for rotational motion compensation and the parameters 
\ hore fm» fo3>---} are the search parameters for rotational motion compensation. The phase 


function of Equation (20) and the rotational angle in Equation (18) gives the relationship 


between rotational angle and dwell time such that [3]: 
O(t)< ©,(t). (21) 
The linear interpolation function in MatLab can reformat the ISAR data so that the radar 


signal is uniformly sampled in @, only linearly dependent on time and the quadratic 


phase error is removed. The final focused ISAR signal, as given in [1], becomes: 





S(t = yA, exp |=) a7 | Ax, + ay,0(0)| (22) 
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where s(x,t'), is the ISAR signal after translational and rotational motion compensation, 
O(t')is the uniformly sampled rotational motion term that only depends linearly on t' and 
t' is the new time variable. 


3. Selection of an Appropriate Basis Function or Phase Function 


In the last section, it was stated that suitable basis functions, h, (t) and h, (¢), 


must be found in order to perform translational and rotational motion compensation. 
There are three methods that can be used to determine if a basis function 1s suitable for 
compensating for motion error. The most common method 1s the projection method 
which it is used in [1], [2], [3] and [6] as the way to determine the suitability of the basis 
functions. The second method follows the work in [16] and [17] by using the STFT and 
the Radon transform to determine basis function suitability by automatically recognizing 
the conditions of focus. The third method of determining basis function suitability uses 
the FFT to automatically recognize focus through the characteristics seen in the power 


spectrum. Note that for these three methods, first the best translational motion compen- 


sation parameters, { f\,,/,>>f;3>---}, are found and translational motion compensation is 


performed on the signal. After this has been completed, the rotational motion compensa- 


tion parameters, { f;,, fy), f53.-.-}, can be found. 
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a. Projection Method to Determine h, (t)and h, (t) 


This approach 1s described mathematically in [1] as: 
hivta his} =argmax|(s(x,)h,, (0) (23) 


where (s(x, t)h,, (1) is the projection of the basis function onto the signal of range bin x 


in the time domain and can be defined, as in [1], as: 
(s(x.h,, ()) =[s(x,2)h%, (sat. (24) 
Rotation motion compensation parameters are found similarly as: 


{Parr faa fos} = arg max |(s(x,0)ph,, (0) (25) 


where (s(x,t)ph - (1))) is defined in a similar fashion as Equation (24). Therefore, from 


the above equations, the most suitable basis function will be the one with the parameters 


\histioe tise} OF {fore Sor» fo3>---} that, after vector multiplication between the signal and 


the basis function and after taking the magnitude of the result, will produce the greatest 
value. Though this method is a very popular method for determining basis function suit- 
ability, it was observed that the projection method did not always lead to an adequate ba- 
sis function. To address this, the next two methods are discussed. 


b. Automatic Recognition of Focus using the STFT and the Radon 
Transform 


The Radon transform has been previously used in ISAR image formation 
in [16] and [17]. In this thesis, the Radon transform was used to determine the suitability 


of the basis and phase functions, , (¢) and ©,(¢), and is a replacement for the previ- 


ously explained projection method. Information on the calculation of the Radon trans- 
form is found in [13] and detection of lines in the time-frequency plane is introduced in 
[5]. This method finds its validity in the argument, as presented in Chapter III, that a fo- 
cused scatterer which has had its motion error and, consequently, its time-varying Dop- 


pler removed, will be a line of zero slope in the time-frequency plane. The first step of 


this method is to take the current search parameters { f,,,/().f;3.---} and formulate the ba- 


23 


sis function, 1, (¢). (If translational motion compensation has already been completed, 
the search parameters { f,,, f5,..f,;,...} are used to form the phase function, ©, (rt), for ro- 
tational motion compensation.) Once /, (¢) or ©, (¢) has been formulated, motion com- 


pensation, either translational or rotational as applicable, 1s performed on the range bin. 

It is only necessary to perform motion compensation on the range bin under examination 
during the basis and phase function searches since the rigid body assumption states that 
motion error is constant from range bin to range bin. The next step is to take the STFT of 


the resulting range bin and to calculate the Radon transform of the resulting image. 


The Radon transform is a proven method for detecting lines in images (see 
[5] and [13]). It is able to return parameterized information on the lines found in the im- 
age. Figure 13 is the plot of the Radon transform of the STFT of range bin 31 in the B727 
simulated data set in its unfocused state, after translational motion compensation and af- 
ter rotational motion compensation. The images of the time-frequency transforms have 
already been displayed for these three states in Figures 11, 12 and 8, respectively. For 
the purposes of this method, only the @ parameter of the Radon transform is important as 
a value of 8 = 90° indicates a flat line in the subject time-frequency plane. Figure 13 


demonstrates that, as lines in the time-frequency plane are flattened to zero slope when 


correct values of {f,,, fis5/,3.--} and {f5,,fs.,fo3,--.} are found and motion compensation 


is correctly applied, the peaks in the Radon transform graphs lie on the 6 = 90° line. The 
final Radon transform, which is taken from the STFT of a focused image shows that there 
are two lines in the image, one laying on the @ = 90° line and the other on the 0 = 89° 


line. Therefore, detecting focus is simply a matter of searching for the parameters 


Whiotio tise} and | f5,,f55,fo3>---} that cause peaks along or very near the 0 = 90° line 


when motion compensation is applied. 


The one significant drawback of this method which unfortunately makes it 
currently unusable, is that both the STFT and the Radon transforms take too long to cal- 
culate and even the most efficient of search algorithms must calculate these transforms 
many times. Therefore, this method cannot be currently used in ISAR image focusing in 


a real-time application. However, with the ever increasing computational power in to- 
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day’s CPUs, it may one day be a suitable for determining basis function suitability. For- 


tunately, the research into this area led to the realization of a method that uses the FFT. 
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Figure 13. Radon Transform of Uncompensated and Compensated Range Bins. 


C. Automatic Recognition of Focus using the FFT 

This method was developed soon after it was realized that the method us- 
ing the STFT and the Radon transform could not be practically implemented. The FFT 
method finds its validity in the explanation given in Chapter III.A in that a time- 
frequency transform of a range bin with motion error will have point scatterers mapped 
out as tilted (or possibly curved) lines and a power spectrum which has areas of maxi- 
mum value smeared across several Doppler bins. A focused range bin will have a scat- 
terers mapped out as lines with zero slope in the time-frequency plane and a power spec- 
trum where maximums are sharp points in there respective Doppler (1.e., cross-range) 


bins (refer to Figures 7 and 8 for images that illustrate this). 


The FFT method starts the same as the previous method in Subsection 


3.B.b of this Chapter, in that the range bin under examination is subjected either to trans- 
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lational or rotational compensation using the test parameters { f,,, fio, f,35---} OF 


{ tre ee Totes , as applicable. The power spectrum of the range bin x 1s then calculated 


using the familiar equations: 


S(f),|_ = FFT(s()r|,)[ (26) 


and 


Sal, = FFT(s(el,)| (27) 


where s(t), ; is the ISAR signal of range bin x after translational motion compensation 
and S(f Yr. is its power spectrum, s(t), i is the ISAR signal of range bin y after rota- 
tional motion compensation and S(/), i is its power spectrum. Now to determine the 


suitability of the basis function, /, (r), or the phase function, ©, (¢), it is necessary to 


determine the following: 


max; S(f),. ; 
thie S ios Aigo-s-f = arg max mie (28) 





Xx 


and 


X40 
bi duisdogsen) = are max max S(S)el,} (29) 


LSP); 


Jy. 


Equations (28) and (29) state that the best search parameters would be ones that place the 
maximum percentage of the power spectrum found in the entire range bin in a single 
Doppler bin. A simple maximum search will not work since incorrect values of the 
search parameters can generate a high peak in a Doppler bin which 1s much greater than 
what will be found if the image were focused, but there will also be a very significant 
amount of the power spectrum in adjacent bins and accordingly the image will be ex- 
tremely blurred. Taking the average is very effective at avoiding this situation because 


the averaging process will score these sets of parameters very low. An example of a 
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range bin, motion compensated with two different sets of parameters, one that does not 
focus the image and the other that focuses the image, is shown in Figure 14. This auto- 
matic recognition of focus technique using the FFT is very fast and is the technique that 
is used in the search algorithms in Chapter IV for determining the suitability of a basis 


function’s or phase function’s search parameters. 
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Figure 14. | Example of a Unfocused and Focused Range Bin and Their Score. 


C. CHAPTER SUMMARY 

In this chapter, two key concepts, which are critical to Chapter IV and the overall 
problem of ISAR focusing, were introduced. First, the mathematical model of the AJTF 
algorithm was presented and it is this algorithm that enables time-varying Doppler due to 
target motion to be removed from the ISAR signal. Secondly, the FFT method of deter- 
mining basis and phase function suitability was presented. Chapter IV uses this method 


as the core of the GA and the PSO search algorithms. 


Zt 
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IV. GENETIC ALGORITHM AND PARTICLE SWARM 
OPTIMIZATION AS SEARCH ALGORITHMS 


This chapter covers the Genetic Algorithm (GA) and the Particle Swarm Optimi- 
zation (PSO) algorithm and their application to extracting translational and rotational mo- 
tion correction parameters from the solution space. The chapter concludes by offering a 
comparison between the two algorithms. 


A. INTRODUCTION TO THE SOLUTION SPACE AND THE EXHAUSTIVE 
SEARCH 


Before beginning the discussion on the GA and the PSO algorithms it is important 
to briefly discuss the construct of the solution space which must be searched and the ex- 


haustive search method, as described in [1] and [3], that 1s typically used to find the 
search parameters { fi,,f 5. / (35: OF {fois fo» Soyo}. In the AJTF algorithm for ISAR 


image focusing, the parameters f,, and f,, are always found using the FFT such that: 


fv DopplerBin jarg max (|FFT (s(2)], Y')} (30) 
and 


f,, = DopplerBin larg max (|FFT (s(t, (31) 


which means that f,, and f,, are always equated to the Doppler bin possessing the max1- 
mum power spectral density in range bins x and y. The remaining parameters, 

{f5,fi3>+-+ and {f5,,.f53,.-.$, must be searched and each parameter can possess any real 
valued number lying between —N to +N, where N is the number of pulses used to form the 
ISAR image in cross-range. Therefore, as the number of pulses that are used to form an 


image increases, so does the solution space and the computational intensity required to 


achieve focus. 


Figures 15 and 16 shows a map of the solution space that must be searched for the 


B727 simulated data set when the search is limited to a two-dimensional search space 


: d ; ‘ : 
(i.e., 3'°-degree motion error which requires a search for parameters { f,,, f,,} and 


| foo» fos). For translational motion compensation, it can be seen that the solution space 
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is fractured, with many areas of separate local maximums. Because of this fracturing, 
conventional searches, such as a gradient search, may have problems finding the global 
maximum and instead converge to a local maximum [9,14]. This local maximum may 
not be an acceptable set of parameters for focusing the image. The rotational motion so- 
lution space is not nearly as fractured but it does have a couple of valleys which may 


cause a gradient search to mistakenly converge to a local maximum. 
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Figure 15. 2D and 3D Visualization of B727 Solution Space for Translational Mo- 
tion Compensation. 
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Solution Space for BY 2? - Rotational Motion Compensation 
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Figure 16. 2D and 3D Visualization of B727 Solution Space for Rotational Motion 
Compensation. 


A common method for searching this solution space, as detailed in [1], [2] and 


[3], is the exhaustive search. The exhaustive search is typically done by finding /,, as in 
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Equation (30) and then stepping through /,, for the span of the solution space (while the 
additional parameters, /,,, f,, and higher terms are set to zero) and finding the best term 


that satisfies Equation (28). Once this term is found, the next term is stepped through the 
solution space, with previous terms set to their optimal values, until again the maximum 
of Equation (28) is found. This is repeated, as required, for all of the search parameters 


until the best basis function, /,, (¢), 1s found and translational motion compensation 1s 


performed on the ISAR signal. This identical methodology is then repeated to find the 
best phase function, ©, (tf). 


Though this is an effective method of finding the parameters of motion compensa- 
tion, it suffers from two drawbacks. The first is that when searching for the third-order 


parameters (i.e., f,, or f,,; and above), the exhaustive search may not find a true global 


maximum since it is doing a 1D search of a parameter with all other parameters fixed, de- 
spite the fact that the parameters are not orthogonal. Secondly, even with this reduced 
search space, it 1s still very computationally intensive and time consuming to step though 
all possible answers. One way of measuring computational intensity is to determine the 
number of times the cost function must be calculated. In terms of this exhaustive search 
algorithm, the cost function is Equation (28) and for the exhaustive search it must be cal- 


culated N 


cost function 


times as determined: 


N 


cost function 


=(2*N 


cross range) NG fisse) NM fpasfiam) (32) 


where N is the number of pulses in the cross-range, N is the number of 


cross _range (fia fiz) 


search parameters for translational motion compensation and N is the number of 


Cae ee 
search parameters for rotational motion compensation. Therefore, for the B727 simulated 
data set with 256 pulses in the cross-range and 3"’-degree motion compensation (required 
to find two search parameters for both translational and rotational motion compensation), 


the number of cost function calculations 1s: 


N 


cost function 


= (2-256)-2-2 = 2048. (33) 


a2 


Therefore, when the efficiency of the GA and PSO algorithm is measured, it will 
be evaluated not only by its run time but also by the number of cost function calculations 
it requires to converge. Minimizing the number of cost function calculations is the key to 
an efficient algorithm and a clear measure of how much better the algorithm is than the 
simple exhaustive search. 

B. GENETIC ALGORITHM PARAMETER SEARCH 

Genetic algorithms, first developed in 1989 by Goldberg, whose work is pub- 
lished in [9], are particularly well suited to searching a large solution space that contains 
one global maximum buried amid many local minima and which may be discontinuous. 
A GA search attempts to model natural behavior where strong traits will propagate 
through a breeding population while weak traits will die off and eventually be removed 
from a population after a number of generations. One outstanding feature of a well- 
written GA is that it will rapidly converge very close to the global maximum. The disad- 
vantage is that once it has converged to this point, typically only small improvements will 
be made to the found global maximum and this only after many generations. This trait is 
illustrated in Figure 17 and is a typical feature of any evolutionary search. The GA ex- 
periences a rapid increase in the best score of its found global maximum in the first 23 
generations of the evolutionary search. The next 32 generations are characterized by sig- 
nificant but smaller increases to the score of the found global maximum. The remaining 
generations are characterized by only a minimal and hardly noticeable increase to the 
score of the found global maximum. Therefore, using this argument, it is obvious to 
state that with GAs, convergence to the absolute and unique maximum 1n not guaranteed. 
However, the probability that convergence has occurred increases with the user’s will- 
ingness to wait for the algorithm to process successive generations. Looking again at 
Figure 17, it can be stated with confidence that after the 70" generation, the GA has con- 


verged either at or near the true global maximum. 


In ISAR imaging, there is no knowledge of the location of the true global maxi- 
mum nor is there any way to predict the likeliness of finding a better global maximum. 
As successive generations occur, the probability of finding a better global maximum di- 
minishes. Also, since the goal is to rapidly find acceptable search parameters that will 
focus the ISAR image, exiting the search as expediently as possible becomes critical. 
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Therefore, the rapid convergence (1.e., increase 1n the score of the found global maxi- 
mum) of the GA early on in the search can be exploited. Combining this rapid conver- 
gence with the ability to optimally decide when to exit the search will allow for rapid fo- 
cusing of ISAR images. The work of Li and Ling in [6] is an excellent first publication 
of a GA applied to ISAR focusing and was used as a reference in the writing of the GA 
code for this thesis. In particular, [6] states that real-valued GAs significantly outperform 
binary coded GAs. Because of this, the GA for this thesis will forego exploring binary- 
coded GAs and investigate only a real-valued GA. Another reference source was [10], 


which provides a good explanation of the flow and construct of a GA. 


Genetic Algorithm 
Best Score and Average Score vs Generation 


Smaller Increases ! 
During Next 32 Generations: 


Very Small Increase 
| /Rapid Increase During After Many Generations 
First 23 Generations 


Score 


—+— Average Score of the Population 
—*— Best Score in the Population 





0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 
Generation 


Figure 17. _ Illustration of a GA’s Rapid Convergence. 


1. Sequence of the Genetic Algorithm 
The genetic algorithm follows a set sequence that repeats until done (see Figure 


18) with each cycle referred to as a generation. The definition of each step 1s as follows: 
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Figure 18. Flow Chart of the ISAR Genetic Algorithm (After [6].) 
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a. Define the Population 

The first step of the GA 1s to generate a population of suitable size to en- 
sure sufficient randomness 1n the initial group of parents. It is important, however, to en- 
sure that the population is not so large that it becomes unwieldy and slows down the algo- 
rithm or that, after only a few generations of the GA, the cost function has been calcu- 
lated more times than it would have been if an exhaustive search had been used. How- 
ever, in the ISAR case, the ideal population size will be function of the number of pulses 
used 1n the cross-range, the number of search parameters due to the order of motion error 
and the maximum number of generations that the GA will run. For the ISAR parameter 
search, the initial population is a matrix with the first element of each parent being f/f, as 
determined by the FFT. (Note that from this point on when search parameters are referred 
to with a single subscribe, the argument applies equally to translational and rotational 
motion compensation search parameters. The double subscript will be added when nec- 
essary to differentiate between the two sets of search parameters.) This parameter does 


not change. The remaining parameters, {/,, f;,...}, are initialized as uniform random 


variables between +N, where N is the number of pulses in the cross-range. After initiali- 


zation, the population will resemble the construct shown in Table 1. 
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Table 1. The Genetic Algorithm Population — 2™ subscript indicated Parent. 


b. Evaluating Fitness of Each Member of Population 

In any search algorithm, the ability to evaluate the suitability of a set of 
parameters compared to other sets of parameters when applied to the cost function is ab- 
solutely crucial in order to ensure that the best parameter set is chosen. In a GA, it 1s 
critical, since the fitness ranking will determine which parents will be used to form the 
next generation and which of the newly formed children will be inserted into the new 
population. Evaluating fitness is the key component in determining the search pattern of 
the solution space. Since this GA uses the FFT version of the automatic recognition of 
focus method to determine parameter fitness, all of the parents or children are evaluated 
against the cost function in Equation (28). In addition, because the calculated fitness be- 
tween superior and poor parameters are small, a scaling exponential (chosen from em- 
pirical observation to be 100) was added to the fitness calculation in order to increase the 
separation between the fitness score of parameters. Care must be taken when choosing 
the exponential. Too low and it will make no difference to the efficiency of the search 
but set 1t too high and the GA will exit part of the solution space and may fail to find the 


global maximum that could reside there. 
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C. Generating Children 

Now that the fitness of each parent has been computed, the next step is to 
form the children that will be part of the next generation. The first step is to determine, 
based on the fitness of each parent, the likeliness that a parent will be mated with another 
parent to form two children. This 1s done using the roulette wheel method as introduced 
in [9] and [10]. The way the roulette wheel works is that each parent is assigned a por- 
tion of the wheel that is proportional to the percentage they contribute to the sum of the 
fitness values of the entire population. This way parents who contribute significantly to 
the population have a significant chance of being chosen to produce children while even a 
poor contributor retains a chance, albeit a small one, to contribute to the population. It is 
critical that even poor parents be given a chance to reproduce since it ensures genetic di- 
versity within the population and leads to a complete search of the solution space. The 
roulette wheel, shown in Figure 19, was built in MatLab by assigning a vector a cumula- 


tive sum of the fitness values (1.e., first value equal 3, second value equal to 13.5 and the 


last equal to 30.02). 
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Figure 19. Illustration of the Roulette Wheel. 
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Once the roulette wheel has been constructed, the mating pool can be con- 
structed. As in Figure 19, if our population is of size 6, we need 6 parents to form the 


mating pool. Choosing is done by generating a pointer which is computed as a linear dis- 


tributed random variable between 0 and LF itness. A parent is chosen if it 1s the first 


parent with a roulette wheel value greater than the random number assigned to the 
pointer. Figure 20 shows a pointer which possesses a value of 19. On the roulette wheel, 
Parent 3 has a value of 18 while Parent 4 has a value of 21. Since 19>18 and 19<21, 
Parent 4 is selected for the mating pool. Parent selection for the mating pool continues in 


this manner until all 6 slots in the mating pool are filled. 


Parent 


Roulette 
Wheel 
Value 





Pointer = 19 In this case, parent 4 is selected 


Figure 20. _ Illustration of the Mating Pool Selection Method. 


Once the mating pool has been populated it is broken into two equal arrays 
and their order 1s randomized. The matching elements of the two breeding pools are the 
mating pairs. Breeding between the two mating pairs are determined by the cross-over 


equations taken from [6]: 
C.=aP,+(1-@B, (34) 
and 


C,=(-a)F,+aF, (35) 
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where F,, and Ff, are the first parents in the two breeding pools, C, and C, are their 


children and q@ is a uniformly distributed random variable between 0 and |. So for each 
search parameter, the child is getting a new parameter that is equivalent to 


Foc = (A) fo ovr +(1-@) i . Also remember that f, does not change as it is de- 


termined by the FFT. 
d. Mutation of the Children 


The purpose of mutation is to allow the search to be global in nature even 
as the parents converge towards each other and a possible global maximum. The muta- 
tion process occurs right after the children have been formed and is illustrated in Figure 


21. Mutations are determined by [6]: 
C, a ac, ae @ — 0) P andom* (36) 


In Equation (36), C, is the child to be mutated, P_,,, 1s a parent pulled randomly from 


the solution space and q@, as before, is a uniform random variable distributed between 0 


and 1. 


Mutation 


Breeding 





Random Parent 1 Child 1 Child 2 Parent 2 
Parent 


Mutated Child 
Solution 
Space 


Figure 21. Illustration of the Mutation Process. 
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The other parameter that needs to be addressed is the mutation rate. Set- 
ting the rate too high destroys the integrity of the local search of a maximum and the GA 
may be terminated before it 1s able to find the absolute maximum. However, setting it 
too low increases the probability that the GA will converge to a local maximum. Refer- 
ence [10] recommends a low mutation rate of 0.1% - 1% since, in nature, mutations are 
relatively rare. Reference [6] uses a mutation rate of 30% for their ISAR GA. Reference 
[11] suggests a variable rate, one which starts high and reduces during the number of cy- 
cles. Empirically, it was determined that the mutation rate from [6], set at a constant 30% 
, produced the best results. In this thesis’ GA, this mutation rate was required to reliably 
break the GA out of a situation where it would repeatedly converged to a local maximum 
that created a poorly focused image. 


é. Evaluating the Fitness of Each Child and Reinsertion of Parents 
and Children into New Population 


Once the mutation process has been completed, each child is evaluated for 
its fitness as described above. The next stage 1s to devise a way to reinsert the old parents 
and the new children into the new population which promotes convergence to the global 
maximum. Reference [10] suggests that reinsertion should be a random process where 
only a certain percentage of the children are randomly chosen to join the new population 
and the remaining empty slots are randomly filled by parents of the old population. Care 
should be taken when applying this approach, particularly when speed of convergence is 
of particular importance. What occurs in this approach is that the current maximum (_.e., 
the best value found so far in the search) could depart the population, poor performers 
will enter into the population and high ranked children may not enter the population. A 
compromise to this, which was initially implemented in this thesis’ GA, was setting the 
reinsertion rate of the children to an 80% probability of insertion into the new population 
with the remaining empty slots in new population being filled by the highest ranked par- 
ents. The probability of convergence graph, as a function of size of population and the 
number of generations evolved using this reinsertion method, was generated and is dis- 
played in Figure 22. (Discussion on how the convergence graph is generated will be dis- 
cussed later in the results section.) As can be seen from this graph, the probability of 
convergence, for 30 trials at each point, has a disappointing maximum value of 80% even 


with a large number of generations and a large population size. To improve this, it 1s 
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tempting to choose a deterministic approach where the best parents and children are cho- 
sen to form the new population. However, while building the GA for this thesis, imple- 
menting reinsertion in this manner tended to lead to convergence to a local maximum and 


not the global maximum. 
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Figure 22. Convergence Graph of GA with Probability of Reinsertion of Children 
into New Population Equal to 80% 


In order to improve the probability of convergence, a simple method 
which combines deterministic and random processes was chosen. When a new popula- 
tion is ready to be formed, 50% of the slots are filled by the best parents of the old popu- 
lation and the best children of the new population. The remaining slots are filled by ran- 
domly choosing from the parents and children not yet chosen for insertion into the new 
population. This method accomplishes two goals. First, it ensures that the best values 


found by the GA will continue to influence the evolution of the GA from generation to 
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generation which encourages local searches around these good parameter sets. Secondly, 
the random selection of the second half of the new population ensures sufficient genetic 
diversity in the population which encourages global searching of the entire solution 
space. This is the method that was employed 1n this thesis’ GA. Another convergence 
graph was generated with this change but it, and its subsequent improvement, are dis- 
cussed in the results section. 

f. Determine if the Search is Complete 

There are a couple of ways to determine if the GA search is complete. 
The first 1s a deterministic way where the user indicates the maximum number of genera- 
tions that the search algorithm is allowed to evolve. After the last generation, the best pa- 
rameter is taken and used for either translational or rotational motion compensation. The 
second method is to measure the number of generations that the best value found has re- 
mained unchanged. The longer this value remains stagnant, the likelihood of finding a 
better value becomes less and less. After a fixed number of generations that this value 
remains stagnant, the GA will stop and use this value for motion compensation. For this 
thesis, the user indicates the number of generations the GA will run since this allows for 
the generation of convergence graphs, which in turn is a good measure of the GA’s per- 


formance and allows the GA to be compared to the performance of the PSO algorithm. 


Finally, this procedure must be completed twice, once for translational 
and then for rotational motion compensation. If compared to the exhaustive search, the 


number of cost function calculations for a GA can be determined by: 


N 


cost function _GA = 2 ; Nice ; Ng (37) 


enerations ° 


Zz: GA ISAR Results 
To determine the functionality of the AJTF GA, the B727 simulated data set and 
several different 2048 pulse imaging intervals of the 6 — Point Delta Wing data set were 
used. 
a. B727 Simulated Data Set 
Shown is Figure 23 1s the unfocused B727 and the power spectras of the 
two range bins used in the AJTF GA. It is easy to see from both the image and the power 


spectras that there 1s time-varying Doppler in the simulated radar data set. Figure 24 
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shows the image after translational motion compensation of 3"-degree motion. Since 
most of the error is rotational, as stated in [3], the image 1s still not focused with the ex- 
ception of one of the two point scatterers in range bin 31, which was changed in its spec- 
trum from being very broad to very narrow. Figure 25 shows the final focused image. 
Notice that the image, after 3""_degree rotational motion compensation, is well focused 
and that the scatterers are now sharp points in the frequency domain. This image was 
formed with a population size of 20 and was evolved for 20 generations for a total run 


time of 0.801 seconds fora N 


cost function _GA 


= 800. As given in Equation (33), an exhaus- 


tive search would have taken 2048 calculations of the cost function. 
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Figure 23. Unfocused B727 and the Power Spectras of Range Bins 31 and 44. 
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Figure 24. B727 and the Power Spectras of Range Bins 31 and 44 after Transla- 
tional Motion Compensation. 
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Figure 25. Focused B727 and the Power Spectras of Range Bins 31 and 44. 


n order to better understand the relationship between the population size, the number of 
generations and the probability of convergence, a convergence graph was calculated. It 
was built by varying the population size, increasing the number of generations that the 
GA would cycle through and, after 50 tests for focus had been made for each point, 
measuring the percentage of times the unfocused image came into focus. Focus is meas- 
ured by taking the 2D cross-correlation of the focused B727 image with the resulting 1m- 
age generated by the GA. Cross-correlating an unfocused image with the focused image 
will create a smeared 2D cross-correlation plane with several points in red as shown in 
Figure 26. Cross-correlating a focused image with the focused image will create a 2D 
cross-correlation plane with one distinct peak as in Figure 27. This can easily be detected 


through threshold detection. 
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Figure 26. Cross-Correlation between a Focused and Unfocused B727 Image. 
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Cross-Carrelation Function of Two Focused Bi?? Images 
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Figure 27. Correlation between Two Focused B727 Images. 


Figure 28 shows the convergence graph as the percentage of images fo- 
cused as a function of varying the population size and the number of generations the al- 
gorithm was allowed to evolve. Figure 29 shows the average run time of each of these 
trials. These graphs are very useful in determining optimal search parameters for the GA 
search. For example, from these graphs it can be seen that an optimal choice of parame- 
ters, for the B727 simulated data set, would be a GA search with population size of 60 
which is evolved for 10 generations since this gives an 84% chance of convergence, a run 
time of only 0.96 seconds and only 1200 calculations of the cost function. The require- 
ment to calculate the cost function 1200 times is only 58.5% of the minimum cost func- 


tion calculations needed by the exhaustive search. 
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Percentage of Trials the Genetic Algorithm Converged to Focused Image 


TOO 
So 
SU 4a 
c=: 
ch 
5 FO% 
ts 
a) 
oe oO 
= 60% 
ay) 
wo 
£ 50% 
O 
> 40% ! 
= 4 ! ! / | —+— 10- Population 
= 30% = : 7 +---| ——+— 16 - Population 
ae | ! ! / | —*— 20 - Population 
20% : ! +---| ——o— 26 - Population 
! ! ( | et 30 - Population 
10% : 7 W--} > 40) - Population 


—P— 60 - Population 





AQ BO 
Number of Generations 


Figure 28. Percentage of Images Focused Using the Genetic Algorithm as a Func- 
tion of Population Size and Number of Generations — B727 Data Set. 


47 
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Figure 29. Average Run Time of the Genetic Algorithm as a Function of Population 
Size and Number of Generations — B727 Data Set. 


b. 6 — Point Delta Wing Experimental Data Set 

The 6 — Point Delta Wing experimental data set, with 2048 (2'') pulses in 
the cross-range, represents a more challenging focusing problem as it does not always 
exhibit ideal motion error characteristics like the B727 simulated data set. Though the 
assumptions of rigid body, two-dimensional motion still, for the most part, hold for this 
data set, there are many imaging intervals for which they do not hold. For the moment, 
the discussion of results will be limited to four imaging intervals and only 2"*-order mo- 


tion error (i.e., a 1-dimensional solution space search for parameter f,, and f,, ) where 


the GA 1s capable of finding a parameter that will focus the image. For comparison, the 
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number of cost function calculations for translational and rotational motion compensation 


using the exhaustive search, as in Equation (32), is 8196 calculations of the cost function. 


The first imaging interval is imaging interval 15 and is shown in Figure 
30. This imaging interval has only one focused point scatterer which can be placed in its 
proper cross-range bin by the fast Fourier transform. The remaining scatterers possess 
time-varying Doppler and are badly blurred in the cross-range. The GA AJTF algorithm 
with a population of size 30 can find the search parameter to focus the image after 50 
generations which equates to 3000 calculations of the cost function or 37% of the calcula- 
tions when compared to the exhaustive search. Notice that the power spectrum in the 
unfocused image 1s spread in Doppler while the focused image has a power spectrum 


with sharp points. 


The second imaging interval is interval 19 and is shown in Figure 31. 
Again, five of the scatterers are out of focus in the cross-range. The AJTF GA is success- 


ful in finding the first order parameters, f,, and f,,, and focusing the image. All six 


point scatterers can be clearly identified in the image. Figure 32 shows imaging interval 
23 whose unfocused image has significant motion error as only one corner reflector is in 
its proper Doppler bin. The focusing algorithm removes the error and focuses the image 
in the cross-range. The sixth reflector is located in range bin 20. However, as it is shad- 
owed by a reflector on the leading edge of the 6 — Point Delta Wing, in range bin 26, it is 
below threshold and not visible. Again both the imaging intervals were formed with 


3000 calculation of the cost function. 


The final imaging interval to be presented is imaging interval 14 in Figure 
33. The original image is very badly blurred in the cross-range. The AJTF GA 1s able to 
focus most of the unfocused scatterers extremely well. There still is some blurring, par- 
ticularly in range bin 25 but also a small amount in range bins 26, 21 and 28. The inabil- 
ity of the AJTF algorithm to remove all blurring is most likely due to a violation of one of 
the basic assumptions used to formulate the mathematical model of the AJTF, most likely 


the presence of a very small degree of 3D motion. 
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Final Image after Rotational Motion Compensation 
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Figure 30. Uncompensated (top) and Motion Compensated 6 — Point Delta Wing 
ISAR Image — Image Interval 15 with 2048 Pulses in the Cross-Range — GA. 
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Final Image after Rotational Motion 
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Figure 31. | Uncompensated (top) and Motion Compensated 6 — Point Delta Wing 
ISAR Image — Imaging Interval 19 with 2048 Pulse in the Cross-Range — GA. 
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Final Image after Rotational Motion Compensation 
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Figure 32. Uncompensated (top) and Motion Compensated Delta Wing ISAR Image 
— Imaging Interval 23 with 2048 Pulses in the Cross-Range — GA. 
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Figure 33. 
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Uncompensated (top) and Motion Compensated Delta Wing ISAR Image 
— Imaging Interval 14 with 2048 Pulses in the Cross-Range — GA. 
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As with the B727 data set, a convergence graph for the Delta Wing data 
set was calculated in order to determine the best combination of population size, number 
of generations in the run time and the resulting probability of convergence. Imaging in- 
terval 15 was used as the test data set to be focused by the GA. The resulting image was 
then compared to the focused version of this imaging interval via the 2D cross-correlation 


function following the methodology as explained for the B727 data set. 


The convergence graph is shown in Figure 34. From the graph, it can be 
seen that for an 80% probability of convergence, a population size of 50 and a run of 50 
generations is necessary. This translates into 5000 calculations of the cost function, or 
61% of what is necessary for the exhaustive search. It is shown in Figure 35 that the av- 
erage run time for a GA with this population size and 50 generational evolutions is 22 


seconds. 
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Figure 34. Percentage of Images Focused Using the Genetic Algorithm as a Func- 
tion of Population Size and Number of Generations: 6 — Point Delta Wing Data 
Set. 


54 


Average Run Time of Genetic Algantim 
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Figure 35. Average Run Time of the Genetic Algorithm as a Function of Population 
Size and Number of Generations: 6 — Point Delta Wing Data Set. 


C. PARTICLE SWARM OPTIMIZATION PARAMETER SEARCH 

The particle swarm optimization (PSO) algorithm was introduced in [7] and finds 
its roots in the fact that the way that fish school and bees swarm is optimal in nature and 
can be adapted for parameter searches to optimally solve engineering problems. In the 


case of the ISAR focusing problem, it is ideally suited to finding the optimal values for 
the search parameters, {f,,, f;;,-..} and {f,,, fi;,...}. as it displays the same search char- 


acteristics as GA, in particular, the ability to rapidly converge towards a global maxi- 
mum. Figure 17 and the arguments that accompany it about initial rapid improvement to 
the found global maximum followed by smaller improvements after many cycles applies 
equally to the PSO as explained for the GA. The driving force behind the particle swarm 
optimization 1s that each particle in the swarm knows the fitness of its current location, 


the best location it has ever encountered, the Local Best, and the best location ever en- 
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countered by any particle of the swarm, the Global Best. It is this knowledge that forms 
the search pattern of the solution space in order to find the optimal solution. 

1. Sequence of the Particle Swarm Optimization Algorithm 

Like the GA, the PSO follows a set sequence that repeats until done. Figure 36 
graphically describes each of the steps in the PSO as it applies to ISAR focusing. The 


following paragraphs describe each step in detail. 
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Figure 36. The PSO Flowchart (After [7].) 


a. Defining the Particles of the Swarm 
The swarm in a PSO 1s identical in definition to the population in a GA, as 
shown in Table 1. Each particle in the swarm 1s analogous to a parent and holds the 


search parameters {/,, f,,/;...}, as determined by the degree of motion error. Each parti- 


cle is initialized as a uniform random variable from +N to —N , where WN is the number 
of pulses in the cross-range. At the same instance of defining the swarm, each particle 1s 
designated an initial velocity that is uniformly random between —2N and +2N. The 


maximum velocity, |v |, that any particle can possess is 2N. Figure 37 shows a 


max 
swarm of a 3"-degree parameter search after initialization. As expected, it is character- 


ized by the uniform spread of the particles through the expanse of the solution space. 
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Initial swarm with Particles 
Randomly Scattered ina 3D Solution Space - 4'_Degree Motion Error 
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Figure 37. Example of an Initial Swarm in a Solution Space for 4""-Degree Motion 
Error. 


b. Assessing fnitial Fitness of the Swarm and Defining the Local 
Best Vector and Global Best Particle 


Now that the initial swarm has been defined, each particle is assessed for 
its fitness as defined by the cost function in Equation (28) 1n a process that is identical to 
that described in the GA section of this chapter. In fact, the MatLab code that defines the 
fitness is identical for the GA and the PSO. Once the fitness is known, the best particle in 
the swarm becomes the Global Best which holds the location of this particle and its fit- 
ness score. The Local Best vector is initialized with the current particle positions and 
their current fitness since this position is each particle’s first position and, therefore, by 
default each particle’s best position they have ever encountered. 

C. Updating the Particle Velocities and Positions 

Now that each of particles has a fitness score assigned to it and its velocity 
initialized, the particle velocities can be updated based on the current location of each in- 
dividual particle, its Local Best particle and the Global Best particle of the swarm. To 
do this Equation (3) in [7] is used: 
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te K(y, TQ, ‘rand *( firgcaiy = I, )+ Q, rand - (Sctobat ay) (38) 


where v, is the particle’s current velocity, f, is the particle’s current location in the ni 
dimension of the solution space and fj, ..4;, ANd foioran are the Local and Global best po- 


sitions, respectively, in the n" dimension. The PSO parameters (~,, ~, and K are pa- 


rameters that, as stated by [7], were determined by extensive research to be optimal when 


set to 2.8, 1.3 and 0.729, respectively. The values of g, and @, controls the balance be- 
tween the attraction of the particle to the local and global bests. With g,> g,, the parti- 


cle is more prone to be attracted towards searching the area surrounding the particle’s lo- 
cal best. The constriction factor, K , can be compared to inertia and it effectively assists 
in the deceleration of the particle as convergence occurs. Rand is a uniformly distributed 
random variable between 0 and | that serves the purpose of adding randomness to this at- 
traction and leads to a more complete spanning of the solution space. We can see from 
Equation (38) that, as the separation between the particle’s current position and the 
swarm’s Global Best and its Local Best expands, the greater effect they will have on the 
particle’s velocity and subsequently its search direction. Also, the search pattern of the 
entire swarm 1s instantaneously affected every time the Global Best changes. Reference 
[7] states that this equation is analogous to a bee’s desire to keep searching for food in the 
direction that it is going but being instinctively pulled towards the best location it has 
ever found (Local Best) and the best location ever found by any member of the swarm 


(Global Best). 


To update the particle’s positions, it is simply a matter of applying Equa- 


tion (39) to each of the search parameters [7]: 
Fn = Sn Vib (39) 


Typically tis taken to be equal to | so, to generate the next cycle of the swarm, it simply 


becomes a matter of adding the velocities, v,, to each appropriate search dimension’s 


current positions. The effect that the velocity update has on the particles is shown in Fig- 


ure 38. 
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Figure 38. Graphical Illustration of the Velocity and Position Update Process. Note 
that gbest = Global Best and pbest = Local Best (From [7].) 


When discussing velocity and position updates, it is necessary to discuss 
the boundary conditions of the solution space since particles left to themselves will depart 
the solution space where meaningful solutions reside. Since we know from previous 


discussion that our search parameters, {/,, f;,...}, will possess values between +JN , each 


dimension of the solution space has a boundary at +N (where WN is again the number of 
pulses in the cross-range). Reference [7] states that there are three ways to deal with par- 
ticles that exceed the boundary of the solution space: 

1. Absorbing Walls: When a particle exceeds a boundary in a dimen- 
sion of the solution space, the velocity in this dimension is immediately zeroed and its 
position is set to the boundary. The particle will then rejoin the swarm when its velocity 
is updated during the next cycle. Extreme care must be taken when using this boundary 
condition since a serious problem that is not mentioned in [7] may occur. If the situation 
exists where a local maximum exists at the boundary and if this position is simultane- 
ously the Local Best and the Global Best position, the velocity will remain stuck at zero 
and the particle will stay fixed to the boundary. A condition was observed where all the 
particles in a swarm became stuck to the boundary and the swarm did not converge to the 
true global maximum and the focused image after motion compensation was very poor. 

ll. Invisible Walls: Here the idea 1s to let the swarm cycle without 1n- 
terruption and any particle that strays beyond the boundaries of the solution space is not 


evaluated against the cost function. Instead it is immediately assigned a fitness equal to 
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zero. As such, this particle outside the boundaries will not encounter any new local or 
global maximums. This will cause the velocity update equation, as in Equation (38), to 
lead the particle back into the solution space without interruption to the particle’s natural 
search pattern. 

111. Rebounding Walls: When a particle passes a boundary, its posi- 
tion is immediately set to the boundary and it is given a velocity that points away from 


the boundary. This new velocity is equal to 


v, =+0.1-rand - Iv 





(40) 


max,, 


where v,_,, 18 the maximum velocity of a particle in the n” direction, rand is the usual 


uniformly distributed random variable between 0 and | and the direction is set such that it 
is heading away from the boundary. The 10% scaling factor keeps the particle near the 
boundary, since that 1s where it was headed, while still allowing the particle to retain 
momentum and avoid the problem of sticking to the boundary, as explained above in the 


rebounding walls boundary condition. 


Reference [7] goes into greater detail on the strengths and weaknesses of 
each boundary condition. However, the PSO algorithm written for this thesis uses the re- 
bounding walls option since it keeps all the particles within the solution space. This cre- 
ates the condition where each particle will have a new position to evaluate against the 
cost function for each cycle, a critical factor when speed of convergence 1s important. 
Also, there is no danger of a problem with the boundary, as explained for the absorbing 


wall boundary condition. 


Finally, when updating the velocity of the particles, it is necessary to en- 


in either the positive or negative direction. Reference 


max 





sure that no particle exceeds Lv 


[7] states that v,,. has been determined to be optimal when set to equal the dynamic 
range of the solution space in a dimension. This is why, as previously stated, 


=2-N. When a particle is observed exceeding the maximum velocity, it is simply 


max 





lv 


set equal to tv_., with the direction of velocity preserved. 


max ? 
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d. Updating Particle Fitness and Reassigning Local Best Vector 
and Global Best Particle 


Now that the new positions of each of the particles have been defined, 
their current fitness values are assessed as previous discussed. If the fitness score of any 
of the particles exceeds that of the Global Best, it becomes the new Global Best. In addi- 
tion, each of the particles current fitness 1s assessed against their individual Local Best. 
If the fitness of their current position exceeds their Local Best, the current position be- 
comes the Local Best. In this way, as new Local Bests and new Global Bests are found, 
the search pattern of each particle and the entire swarm is updated to bias their search to- 
wards these areas. 

é. Termination of the Search 

At the termination of the search, the swarm should resemble the swarm 
shown in Figure 39 and should have converged to the solution space’s global maximum. 
The PSO is terminated after a fixed number of cycles as chosen by the user, which is the 
same way that this thesis’ GA is terminated. As with the GA, terminating the search in 
this manner adds a deterministic quantity to a random search which allows the PSO to be 
evaluated for performance and compared to the GA 1n terms of probability of conver- 
gence. Similar to this thesis’ GA, on completion of the last cycle, the Global Best parti- 
cle is taken and used for translational or rotational motion compensation as applicable. 
Also, the number of cost function calculations of the PSO algorithm can be determined 
using Equation (37), where the variable names have been altered to conform to the con- 


ventions of a PSO: 


N = ee eee (41) 


cost function _ Swarm particles 
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Final Swarm with Particles 
Randomly Scattered in a 3D Solution Space - 4'"_Degree Motion Error 
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Figure 39. Example of a Final Swarm in a Solution Space for 4'"-Degree Motion 
Error. 


ye PSO ISAR Results 
To determine the functionality of the AJTF PSO algorithm, the B727 simulated 
data set and the four 2048 pulse imaging intervals of the 6 — Point Delta Wing experi- 
mental data set, were used. These four imaging intervals were the same intervals that 
were used in the GA ISAR results. 
a. B727 Simulated Data Set 
A focused B727 image is found in Figure 40 with the spectra of the range 
bins used for focusing shown on the right of the image. It was focused by the PSO algo- 
rithm in 0.681 seconds after 20 cycles and with a swarm size of 20, or 800 calculations of 
the cost function. As with the GA, it is much more important how this algorithm behaves 
in a consecutive series of runs rather than just in one run instance. Using the PSO algo- 
rithm, a convergence and run time graph was constructed to verify algorithm perform- 
ance when tasked to focus this simulated data set. The convergence graph is shown in 
Figure 41 while the average run time graph is shown in Figure 42. These two graphs 
show that an optimal selection of parameters could be 10 particles and 50 cycles as this 
produces an 88% chance of convergence, a run time that 1s below | second and requires 
49% of the required cost function calculations of an exhaustive search. If probability of 
convergence was of critical importance, the parameters of 26 particles and 30 cycles 


could be chosen. These parameters yield a 94% chance of convergence at a cost of a run 
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time of 1.2 seconds and 76.1% of the cost function calculations as compared to the ex- 


haustive search. 
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B727 Data Set Focused by the PSO Algorithm. 
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Figure 41. Percentage of Images Focused Using the Particle Swarm Optimization 
Algorithm as a Function of Swarm Size and Number of Cycles — B727 Data Set. 
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Average Run Time of the Particle Swarm Algorithm 
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Figure 42. Average Run Time of the Particle Swarm Optimization Algorithm as a 
Function of Swarm Size and Number of Cycles. 


b. 6 — Point Delta Wing Experimental Data Set 

The Particle Swarm Algorithm was tested against four 2048-pulse imaging 
intervals of the 6 — Point Delta Wing experimental data set in order to ensure that it could 
focus real radar data. These four imaging intervals, 15, 19, 23 and 14, are the same inter- 
vals that were used in the GA section of this chapter and are shown in Figures 43 — 46. 
The focused images, which were generated by the PSO algorithm, are exactly the same as 
those generated by the GA. The images were typically formed with 1200 calculations of 
the cost function over a run time of 6 seconds. Again, to better understand the relation- 
ship between swarm size, number of particles and their effect on the probability of con- 
vergence and run time, a convergence graph and a run time graph were constructed and 


are displayed in Figures 47 and 48. 
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Final Image after Rotational Motion Compensation 


Delta Ving Data Set - Imaging Interval 15 - Power Spectrum 
¥ 10 Range Bin 25 









18 10 
3 
hh 
TW 
20 = 6 
= 
= 4 
= 
ro 2 
a 
= a 
o 1010 101s: - 20 1025- 1030> 1035 
i A 
: Cross Range - Doppler Bins 
= Power Spectrum 
fe x ‘Tai Range Bin 20 
= 20 
= 
_ 
Oo 10 
20 oo & 
a 
a 
a B 
30) rt 4 
eC 
32 —_ — = . 
(OU TOs T0207 025. 10s0- 10as 1010 1015 1020 10275 1030 1035 
Cross Range - Doppler Bins Cross Range - Doppler Bins 


Figure 43. Motion Compensated 6 — Point Delta Wing ISAR Image — Image Inter- 
val 15 with 2048 Pulses in the Cross-Range — PSO. 
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Figure 44. Motion Compensated 6 — Point Delta Wing ISAR Image — Image Inter- 
val 19 with 2048 Pulses in the Cross-Range — PSO. 
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Final Image after Rotational Motion Compensation 
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Figure 45. Motion Compensated 6 — Point Delta Wing ISAR Image — Image Inter- 
val 23 with 2048 Pulses in the Cross-Range — PSO. 
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Figure 46. Motion Compensated 6 — Point Delta Wing ISAR Image — Image Inter- 
val 14 with 2048 Pulses in the Cross-Range — PSO. 
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From Figures 47 and 48, it can be seen that the PSO algorithm performs 
very well against the 6 — Point Delta Wing experimental data set and for an 84% prob- 
ability of convergence, it requires an average run time of only 14.3 seconds and 3000 cal- 
culations of the cost function. This represents 37% of the cost function calculations re- 
quired by the exhaustive search. Also, Figures 47 and 48 shows that the PSO is capable 
of performing extremely well, even with a small number of particles in the swarm. For 
example, using 20 particles and running the algorithm for 50 cycles leads to a total of 
2000 calculations of the cost function (24% of the exhaustive search), a 9-second run 
time and a 98% chance of convergence to focus. This shows that this algorithm pos- 
sesses both speed and reliability even against a real radar image with a large number of 


pulses in the cross-range. 
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Figure 47. Percentage of Images Focused Using the Particle Swarm Algorithm as a 
Function of Swarm Size and Number of Cycles — Delta Wing Data Set. 
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Average Run Time of Particle Swarm Algorithm 
Delta Ving Data Set - 2" Order Motion - 50 Trials 
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Figure 48. Average Run Time of the Particle Swarm Optimization Algorithm as a 
Function of Swarm Size and Number of Cycles — Delta Wing Data Set. 


D. COMPARISON BETWEEN GA AND PSO ISAR FOCUSING 

ALGORITHMS 

Comparisons between the GA and the PSO algorithms are made in [7] and [11] 
and they typically find that both evolutionary search techniques produce similar results. 
In particular, [11] states that during some optimization trials, which in this case was op- 
timizing the complex weights of an antenna array, the GA outperforms the PSO while in 
others the PSO outperforms the GA. In the case of the ISAR focusing problem, for imag- 
ing intervals where the AJTF algorithm was a valid mathematical model, the resulting fo- 
cused images, produced by the GA and PSO algorithms, were of equal quality as both 
search techniques were capable of finding the optimal search parameters. When a poorly 
focused image was produced by the GA or PSO, it was from an imaging interval of the 6 
— Point Delta Wing experimental data set that possessed qualities that were in severe vio- 
lation of the assumptions of the AJTF motion compensation model. These imaging inter- 
vals were of equivalent poor quality regardless of whether the GA or the PSO was used 
for the focusing parameter search. Also, by examining the average run time graphs of 
both the GA and the PSO, it can be concluded that, for equivalent calculations of the cost 
function, the two algorithms are equally fast at processing through the user-defined num- 


ber of generations or cycles. However, if the convergence graphs for the B727 simulated 
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data set in Figures 28 and 41 and the convergence graphs for the 6 — Point Delta Wing 
experimental data set in Figures 34 and 47 are examined, it 1s shown that there is a sig- 
nificant difference in performance between the two search algorithms. Figure 49 com- 
bines part of the GA and the PSO results from their respective convergence graphs with 
the lines with triangles denoting GA results and lines with circles denoting PSO results. 
What can be seen from the graphs in Figure 49 is that the PSO written for this thesis con- 
sistently outperforms the GA by reliably converging to a focused image even when the 
PSO has less particles in the swarm than the GA has parents in its population. For the 
B727 data set, a PSO with 10 particles outperforms a GA with a population size of 10 and 
16 and performs almost as well as a GA with a population size of 26. In the 6 — Point 
Delta Wing experimental data set, a PSO with 20 particles outperforms a GA with a 
population size of 20, 30 and 50. Also, from a programming point of view, a PSO is 
much easier to code and to troubleshoot and its search pattern is judged merely by Equa- 
tion (38) and, although some of the constants given by [7] can be changed, there was no 
requirement to implement any parameter changes for the PSO in this thesis. For a GA, 
setting a correct mutation rate, reinsertion rate and selecting parents correctly for breed- 
ing are all of critical importance for reliable convergence and can take considerable effort 
to set at a correct level that balances global and local searches. For these reasons, the 
PSO is determined to be better suited for ISAR focusing and only the PSO will be used 


for the subsequent sections of this thesis. 
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Figure 49. Combined GA (top) and PSO Convergence Graphs for Selected Popula- 
tion / Swarm Sizes. 
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E. HIGHER ORDER MOTION COMPENSATION 


Higher order motion occurs when a target’s motion 1s extremely chaotic. This 


creates the situation where searching for the 2"°-order search parameters, f,, and f,,, is 


insufficient to focus the ISAR image. This has been partially addressed already in the 


B727 simulated data set where it was required to find parameters {f,,, f,,} and 
{ Toa fy;} to focus the image. However, with only 256 pulses in the cross-range, which 


translates to an ideal solution space of 262,144 possible combinations, the PSO 1s able to 
quickly traverse the solution space and find the focusing parameters. The problem with 
higher-order motion compensation occurs when there are a large number of pulses in the 
cross-range. For example, for 3"-order motion compensation, the solution space of the 


2048 pulse 6 — Point Delta Wing experimental data set increases from 4096 possibilities 
to 16.7x10° possibilities. If 4"_order motion compensation is required, the solution 


space expands to 68.710” possibilities. Also, do not forget, that these solution spaces 
must always be traversed twice, once to extract translational and the other to extract rota- 
tional motion compensation parameters. This creates an enormous search space that can- 


not be quickly traversed by even the most efficient PSO. 


Imaging interval 13 of the 6 — Point Delta Wing experimental data set is given as 
an example and its unfocused image is shown in Figure 50. This image possesses 4""- 
order motion error and cannot be focused with 2""-order motion compensation. The exis- 
tence of higher-order motion can be seen by examining a range bin in the time-frequency 
plane. As explained in [3], scatterers with higher-order motion will be curved as com- 
pared to 2"’-order motion error whose scatterers are tilted in the time-frequency plane. 
The time-frequency transforms of the scatterers in the unfocused range bins 25 and 19 are 
shown on the left of Figure 51. The fact that they are curved is a strong indication of the 
presence of higher-order motion error. If the standard PSO AJTF focusing algorithm, us- 
ing 4"_order motion compensation, is applied to the image, an extremely well-focused 
image is produced as seen in Figure 52. Also, the right side of Figure 51 shows that the 
time-frequency transforms of the scatterers from the focused image have been straight- 


ened since the time-varying Doppler has been removed from the signal. Unfortunately, 


this requires approximately 110° calculations of the cost function. This limits the PSO 
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usefulness since this requires an excessive run time of | hour and 10 minutes. The way 
around this is to adapt the AJTF PSO algorithm to perform the search in the same pattern 
as the exhaustive search, which was introduced earlier in this Chapter. In this variation 


and, like in the exhaustive search, the PSO is applied to searching for f,. Once the op- 
timal fit is found, /f, is searched for with the PSO. Once the optimal /, is found, as long 
as the addition of the /, parameter adds to the fitness value as calculated by Equation 
(28), it is added to the basis function. If the fitness decreases, f, 1s set equal to zero. 


This process is continued for as many degrees of motion error as required and is com- 
pleted for both translational and rotational motion compensation. The resulting focused 
image, for 2"¢ 3" and 4"-order motion compensation is shown in Figure 53. The ae 
order motion compensated image is shown on the left and though it is not as good as the 
image in Figure 52, it can be achieved in 20 seconds with only 4500 calculations of the 
cost function. Therefore, a sub-optimal, but usable, image can be achieved from an im- 
age blurred by higher-order motion using this focusing method which is exceptionally 


fast even with a large number of pulses in the cross-range. 
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Figure 50. Unfocused Imaging Interval 13. 
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Figure 51. |Time-Frequency Transforms of Range Bins 25 and 19 — Left: Unfocused 
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Figure 52. Focused Imaging Interval 13. 
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Figure 53. Image Interval 13 with 2", 3“ and 4"-Degree Motion Compensation, 
Reduced Solution Space. 


F. CHAPTER SUMMARY 

In this chapter, we have examined the construct of the solution space in the ISAR 
focusing problem and saw that it possessed many local maxima but only one global maxi- 
mum. This characteristic makes the solution space ill-suited for conventional searches 
but well suited for the GA and the PSO evolutionary searches. It was also shown that the 
GA and the PSO provided considerable computational time saving over the exhaustive 
search which is the method normally used for parameter extraction in the AJTF algo- 
rithm. Also, it was shown that the PSO performed considerably better than the GA in the 
ISAR focusing task. Finally, higher-order motion compensation, using the PSO AJFT, 
was demonstrated on an imaging interval of the 6 — Point Delta Wing experimental data 
set. In the next chapter, we will examine what occurs when 3D motion enters the radar 
signal and how this can be accommodated since this is a violation of a fundamental as- 


sumption of the AJTF. 
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Vv. 3D MOTION DETECTION 


In Chapter IV, a GA and PSO algorithm were successfully adapted to the optimi- 
zation of the AJTF algorithm which provided a solution to the computational intensity 
required for successful ISAR focusing that was one of the algorithms major weaknesses. 
A number of imaging intervals, which are now well focused in the cross-range, were pre- 
sented. These imaging intervals, however, follow the mathematical model of the AJTF. 
The next step in this thesis is to examine an effective method for dealing with imaging in- 
tervals which violate the AJTF’s mathematical model. 


A. INTRODUCTION TO THE 3D MOTION MODEL AND THE 
DETECTION METHOD 


One of the basic assumptions of the adaptive joint time-frequency motion com- 
pensation model is that the rotational motion used to form the ISAR image is confined to 


a 2D plane during the coherent processing interval, T.,,.. When motion in the roll plane 


of the aircraft occurs during the coherent processing interval, as shown in Figure 54, 3D 
motion is present in the ISAR return signal and the motion error cannot be compensated 
by the motion compensation model in its current form. Due to the difficulty in compen- 
sating for 3D motion error, Chen, Li and Ling, 1n [3] and [4], developed the linearity of 
phases method to detect the presence of 3D motion. Thayaparan further expanded on the 
3D detection method in [12] and it 1s their work that was adapted to the AJTF PSO mo- 
tion compensation algorithm for the purposes of 3D motion detection in this thesis. The 
goal of 3D motion detection is to determine which imaging intervals of an ISAR data set 
contains a low degree of 3D motion and can be focused using the AJTF algorithm. Imag- 
ing intervals found to possess a high degree of 3D motion are simply discarded and not 


presented to the ISAR operator. 
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Figure 54. Engagement of a Target with 3D Motion (From [3].) 


Now that there is 3D motion in the ISAR signal, Equation (1) can be written in a 


more general form of [3]: 
= ey, 
s(x,t)=) A, exp] =i [RO + + HAO + 2.00] (42) 


where @(t¢) describes the independent angular parameter in the roll plane and Cm Vee.) 


are the 3D coordinate position of the k” scatterer. Chen explains in [3] that Equation (42) 
will reduce to the 2D model in two cases. The first case is when z, is small and the ¢(¢) 
phase term can be neglected. The second case is when there exists a linear relationship 


between @(¢) and @(¢) such that @(¢)=a0(t) [3]. 


The first step towards detecting 3D motion is to remove translation motion error 


from the signal by finding the appropriate basis function, /, (t), as previously described 


and multiplying the signal by its conjugate. With the translational motion error removed, 


Equation (42) reduces to [4]: 
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ee t), = =A exp |=) “hoax ax, +4y,0(0)+85,0(0]] 
(43) 


-YA exp] La —-°| Ax + (Ay, + aAz ao 


Now that the translational motion error has been compensated, 1t becomes necessary to 
find the rotational motion compensation search parameters, { f),, f5.,/53.-..}, that are ap- 
propriate to formulate the phase function, ©, (¢), as in Equation (20), for several differ- 


ent point scatterers in the target. For the purposes of this thesis, 3 point scatterers from 
the 6 — point Delta Wing data set are used and the AJTF PSO that was developed above is 


adapted to extract the search parameters. 


Once all three phase functions have been extracted we have: 


| l 
On) aJilt aioe taal hes 
2! 3! 
It | 
Op. (t) = fa,t+ 5, Saat +o fat ta, (44) 
It 
©,3(t) = fat 5 Sat aga +... 


For signals that do not possess 3D motion, any of the above phase functions, ©,, (t), 
©,,(t¢) or ©,;(¢), would be suitable for rotational motion compensation which would 


create an ISAR image that is well focused in the cross-range. For 3D motion detection, 


Thayaparan’s work in [12] is followed and an average phase function is formed: 


Fo, + far, + far, eerie + fro, + foo, Je +o {& ES eEEN wee Je ee (45) 
2! 2 


The purpose to formulating © (¢) is that as more and more search parameters are ex- 


avg 
tracted from more and more point scatterers suitable for rotational motion compensation, 


Oi (¢) vere’ + @. (4) where ©,,,,,(¢) is the best phase function possible that de- 


scribes the rotational motion error in the whole signal and, after rotational motion com- 


pensation, produces the best image. In this discussion, the number of point scatterers 1s 
Va 


limited to three since the maximum number of scatterers in the 6 — Point Delta Wing data 


set is six and not all scatterers are suitable for rotational motion compensation. 


Now, if the imaging interval possesses a low degree of non-linearity, the plot of 


©,,.0,, or ©,, as a function of O,,,, would be a straight line as in Figure 55. There- 


fore, the first step in determining the degree of non-linearity 1s to find the straight line of 
best fit, in a least-squares sense using the MatLab polyfit function, between ©,,,0,, or 
©,, and O,,,,,. The second step is to plot ©,,,0,, or ©,, as a function of ©.,,,,, and 


measure the percentage that this plot deviates from the line of least-squares fit using: 


_ . Op: (Gece (t)) — LSF{®p, (jeu (0))} 
mt SFO (ies ()P ” 


where NF, is the degree of non-linearity between ©,,(t) and ©... (t), Ox: (Qsdea (t)) is 


the plot of ©,,(¢) as a function of O,,.., (4), LSF{O,, (© sical (r))} is the straight line 

least-squares fit between ©,,(¢) and O,,,,, (4) and N is the number of pulses in the cross- 
range. This procedure is repeated to generate NF, and NF; for the other two phase func- 
tions and the imaging interval is assigned its degree of non-linearity by simply averaging 


the non-linearity functions [12]: 


NF, + NF, + NF- 
ONE preec _ —._ ° (47) 
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Plot of ©,, as a Function of ©,,,,, when no 3D motion is present (After 
[12].) 
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Figure 55. 


Figure 56 shows the Non-Linearity of Phases plot for an imaging interval with a 
high degree of non-linearity. Note that the plot of Op; (Oe (t)) deviates significantly 
from the line of least-squares fit. Figure 57 shows the Non-Linearity of Phases plot for 
an imaging interval with a very low degree of non-linearity. Note that the plot of 


Op: (Oj4eai(t)) in this case is plotted exactly on top of the line of least-squares fit. 
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Figure 56. Imaging Interval Possessing a High Degree of Non-Linearity (After 
[12].) 
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Figure 57. Imaging Interval Possessing a Low Degree of Non-Linearity (After [12].) 
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B. 3D MOTION DETECTION RESULTS 

The 3D motion detection algorithm was tested against the 6 — Point Delta Wing 
data set with 2048 and 1024 pulses in the cross-range and 2” and 4"-degree motion com- 
pensation. The primary goal of this section is to demonstrate the capability of Thaya- 
paran’s 3D motion detection algorithm to differentiate between imaging intervals pos- 
sessing a high degree of non-linearity and those that possess a low degree of non-linearity 
after the algorithm has been adapted to the AJTF PSO algorithm,. Secondly, the effect of 
choosing different sizes of imaging intervals and different degrees of motion compensa- 
tion in detecting 3D motion will be examined. 


1. Imaging Intervals with 2048 Pulses and 2"°-Degree Motion Compen- 
sation 


If the 3D motion detection algorithm is applied to the 6 — Point Delta Wing ex- 
perimental data set with its imaging intervals divided into 2048 pulses and 2""-degree 


motion compensation 1s applied, the non-linearity of phases graph 1s shown in Figure 58. 
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Figure 58. Degree of Non-Linearity of Imaging Intervals of 2048 Pulses and 
2"™'-Degree Motion Compensation. 
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Figure 58 indicates that imaging intervals with a low degree of non-linearity of 
phases are imaging intervals 23, 15, and 18. If these images are focused using the AJTF 


PSO designed for this thesis, the resulting well focused images are shown in Figures 59 — 
61: 
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Figure 59. Imaging Interval 23 — Well Focused with Low Degree of Non-Linearity 
of Phases. 
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Figure 60. Imaging Interval 15 — Well Focused with Low Degree of Non-Linearity 
of Phases. 
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Figure 61. Imaging Interval 18 — Well Focused with Low Degree of Non-Linearity 
of Phases. 
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Figure 58 also indicates that imaging intervals 11, 12 and 26 possess a high de- 
gree of non-linearity of phases. If these images are focused using the AJTF PSO de- 
signed for this thesis, the following poorly focused images, shown in Figures 62 - 64 re- 


sult: 
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Figure 62. Imaging Interval 11 — Poorly Focused with High Degree of Non- 
Linearity of Phases. 
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Figure 63. Imaging Interval 12 — Poorly Focused with High Degree of Non- 
Linearity of Phases. 
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Figure 64. Imaging Interval 26 — Poorly Focused with High Degree of Non- 
Linearity of Phases. 
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Zs Imaging Intervals with 2048 Pulses and 4""-Degree Motion Compensa- 
tion 


Now, if 4"_degree motion compensation (1.e., motion compensation parameters 
\ hove So» So3s fog} ) are used, the following non-linearity of phases graph is generated and 


shown in Figure 65: 
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Figure 65. Degree of Non-Linearity of Imaging Intervals of 2048 Pulses and 4"- 
Degree Motion Compensation. 

In general, there 1s consistency in that those imaging intervals which had a low 
degree of non-linearity of phases in Figure 58, still register as having a low degree of 
non-linearity in Figure 65. The same applies for the imaging intervals with a high degree 
of non-linearity. There are two reasons for the differences between Figures 58 and 65. 
First, the AJTF PSO used to extract the search parameters has a certain randomness asso- 
ciated with its search pattern and it is possible that, for certain imaging intervals, incor- 


rect search parameters were extracted and thus the degree of non-linearity calculation was 
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slightly off. Secondly, 4""-degree motion compensation is more accurate so, provided 
that all of the search parameters were extracted correctly, 1t provides a better measure on 


the degree of non-linearity. 


Figure 65 indicates that the best three imaging intervals are imaging intervals 15, 
14 and 23. Imaging intervals 15 and 23 are the same as before and are shown above. 
Imaging interval 14 is new and 1s now determined to be one of the best three imaging in- 
tervals in the data set. If this image in Figure 66 1s compared to imaging interval 18 in 
Figure 61, it can be seen that imaging interval 14 1s a better image than the one produced 


by imaging interval 18. 
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Figure 66. Imaging Interval 14 — Well Focused with Low Degree of Non-Linearity 
of Phases. 


Figure 65 indicates that the three worst imaging intervals are imaging intervals 
26, 11 and 28. Imaging interval 12, which 1s a particularly bad imaging interval and has 


been shown before in Figure 63, falls off the list of the top 3 worst imaging intervals. 
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Fortunately, it is still ranked as having a high degree of non-linearity and would still not 


be imaged. The new imaging interval, 28, is shown in Figure 67 and is not very well fo- 


cused in the cross-range. 
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Figure 67. Imaging Interval 28 — Poorly Focused with High Degree of Non- 


ae 


Linearity of Phases. 


Imaging Intervals with 1024 Pulses and 2"°-Degree Motion Compen- 
sation 


There are several reasons to use imaging intervals with fewer pulses in the cross- 


range. The most obvious is that is that the 7,,, 1s smaller and subsequently there should 


be less motion error in each imaging interval. Secondly, since there are more imaging 1n- 


tervals available for processing, if an imaging interval has to be discarded, there are still 


many more imaging intervals available to work with. The downside to this is that with 


less pulses in the cross-range, the cross-range resolution, Ar, 


will be less. 


ri? 
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Figure 68 shows the non-linearity of phase graph for 1024 pulses in the cross- 
range and 2"’-degree motion compensation. The three best imaging intervals are imaging 
interval 45, 33 and 35 and are shown below in Figures 69 — 71. The three worst imaging 
intervals are imaging intervals 56, 25 and 29 and are shown 1n Figures 72 — 74. Again, 
imaging intervals with a low degree of non-linearity are easily focused while imaging in- 


tervals with high degree of non-linearity do not focus. 
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Figure 68. Degree of Non-Linearity of Imaging Intervals of 1024 Pulses and 2""- 
Degree Motion Compensation. 
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Figure 69. Imaging Interval 45 — Well Focused with Low Degree of Non-Linearity 
of Phases. 
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Figure 70. Imaging Interval 33 — Well Focused with Low Degree of Non-Linearity 
of Phases. 
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Imaging Interval 35 Low Non - Linearity 
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Figure 71. Imaging Interval 35 — Well Focused with Low Degree of Non-Linearity 
of Phases. 
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Figure 72. Imaging Interval 56 — Poorly Focused with High Degree of Non- 
Linearity of Phases. 
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Imaging Interval 25 High Non - Linearity 
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Figure 73. Imaging Interval 25 — Poorly Focused with High Degree of Non- 
Linearity of Phases. 
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Figure 74. Imaging Interval 29 — Poorly Focused with High Degree of Non- 
Linearity of Phases. 


9] 


4. Imaging Intervals with 1024 Pulses and 4""-Degree Motion Compensa- 
tion 


If 3D motion detection is repeated for 1024 pulse imaging intervals but this time 
using 4""-degree motion compensation, the degree on non-linearity graph is generated as 
in Figure 75. It can be seen from this graph that imaging intervals 10, 33 and 2 possess a 
low degree of non-linearity while imaging intervals 24, 56 and 29 possess a high degree 
of non-linearity. Good imaging intervals 10 and 2 are shown below in Figures 76 and 77 
and poor imaging interval 24 is shown in Figure 78. Again, there 1s consistency between 
the 2" and 4""-degree linearity of phase graphs (in Figures 68 and 75). Imaging intervals 
that register as having a high degree of non-linearity of phases for the 2"-degree motion 
compensation also register as having a high degree of non-linearity of phases for 4""- 


degree motion compensation. 


Degree of Non-Lineanty in Each Imaging Interval 
6 - Point Delta Wving Data Set with 1024 Pulses and Degree of Motion Compensation - 4 


800 
FOO 


600 


Degree of Nan-Linearity 


200 


100 





U 10 Zt) a) At) OU 
Imaging Interval 


Figure 75. Degree of Non-Linearity of Imaging Intervals of 1024 Pulses and 4"- 
Degree Motion Compensation. 
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Imaging Interval 10 Low Non - Linearity 
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Figure 76. Imaging Interval 10 — Well Focused with Low Degree of Non-Linearity 
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Figure 77. Imaging Interval 2 — Well Focused with Low Degree of Non-Linearity of 
Phases. 
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Figure 78. Imaging Interval 24 — Poorly Focused with High Degree of Non- 
Linearity of Phases. 


5. Comparison between the Different Imaging Interval Sizes and the De- 
gree of Motion Compensation 

Now that the degree of linearity of phases has been generated for the two different 
imaging interval sizes and 2" and 4""-degree motion compensation, it is possible to com- 
pare the results and infer which combination is better suited for rapid ISAR image focus- 
ing. The most straight-forward way to go about this is in the form of a table that lists the 
processing time required to produce each graph. Note that the total processing time is the 
time that it takes to bring every imaging interval in the 60,000 pulse 6 — Point Delta Wing 
experimental data set from their non-focused state to the point were their degree of 3D 
motion has been measured and a decision can be made with respect to whether or not the 


imaging interval should be focused or discarded. The results are shown in Table 2: 
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# of Pulses in Total Process- Processing 


ing Time Time per Image 
Cross-range 


(seconds) (seconds) 





Table 2. Processing Time — Degree of Non-Linearity of Phases for Different Im- 
aging Intervals. 


The first thing that is noticeable from Table 2 is that, unsurprisingly, the process- 
ing time for the entire data set increases by over a factor of 3 when 4""-degree motion 
compensation 1s used instead of 2"'_degree motion compensation. Since there is only a 
marginal increase in the quality of the images (only imaging interval 13 with 2048 pulses 
is significantly affected), using higher order motion compensation for the 6 — Point Delta 


Wing experimental data set is a poor use of processing time. 


This now leaves the 2048 pulse imaging interval with 2"-degree motion compen- 
sation to be compared its 1024 pulse counterpart. Using the data set divided into 1024 
pulse imaging intervals has two advantages over the 2048 pulse imaging intervals. First, 
it has twice as many images to work with so, when the AJTF PSO does not converge to a 
focused image (which will occasionally be the case when the PSO 1s set at an optimal 
combination of probability of convergence and run time), discarding an imaging interval 
is no great loss. Secondly, half of the processing time per image is helpful in getting im- 
ages to the ISAR operator as quickly as possible. The data set with 2048 pulses in the 
cross-range has the advantage of producing better images with point scatterers that are 
well separated in Doppler. In an operational situation, these trade-offs will have be bal- 
anced against each other. However, for the purposes of this thesis 1t was easier to set up 
3D motion detection for the 2048 imaging intervals since 3D motion detection and the 
AJTF algorithm are very sensitive not only to range bin selection but as well to point 


scatterer selection. 
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C: CHAPTER SUMMARY 

This Chapter presented Thayaparan’s method of 3D motion detection after it had 
been adapted to the AJTF PSO. It was very effective and efficient at sorting between the 
good imaging interval which focused well and the poor imaging intervals that possessed 
3D motion and could not be focused. The poor imaging intervals can simply be dis- 
carded and not presented to the ISAR operator. In addition, 1t was shown that for the 6 — 
Point Delta Wing experimental data set, 2""-degree motion compensation is the most effi- 


cient used of computational time. 
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VI. CONCLUSION 


A. SUMMARY OF RESULTS 

1. AJTF GA and PSO Algorithms 

The AJTF GA and PSO algorithms performed very well in their task of reducing 
the computational burden required by the AJTF algorithm when employing the exhaus- 
tive search for search parameter extraction. Typically, less than 50% of the cost function 
calculations were required for the GA and PSO algorithms to extract the translational and 
rotational motion compensation parameters which were able to focus the ISAR images in 
the cross-range. In addition, a new method of determining search parameter suitability, 
the automatic recognition of focus using the fast Fourier transform, was designed and 
proved very fast and extremely effective in providing accurate cost function returns that 
gave an accurate figure of merit on the suitability of solution space search parameters. 
When comparing the two different algorithms, the PSO algorithm consistently outper- 
formed the GA and was able to converge to a focused image with less calculations of the 
cost function. Also of important note 1s that the use of the PSO algorithm in ISAR image 
focusing is a unique application of this evolutionary search technique which has not been 
published in the scientific literature. 

Zs 3D Motion Detection 

Thayaparan’s 3D motion detection algorithm in [12], after adaptation to this the- 
sis’ AJTF PSO algorithm and uniquely changing the measure of non-linearity of phases 
to the percentage of deviation from the line of least-squares fit, proved effective in sepa- 
rating good imaging intervals, possessing only 2D motion, from the poor imaging inter- 
vals that possessed a significant amount of 3D motion. Application of the AJTF PSO to 
the 3D detection algorithm significantly reduced computational time, especially since it 
was required that 3 phase functions be extracted from 3 different point scatterers. The 
AJTF PSO algorithm is a unique approach to the 3D motion detection problem and cre- 
ates an algorithm which 1s both fast and accurate. It was also realized that the applica- 
tion of motion compensation greater than the 2""-degree just adds to the computational 
burden associated with 3D motion detection without adding any significant improvement 


to 3D motion detection or image quality. 
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B. FOLLOW-ON WORK 

For those interested in ISAR image formation, there are several areas of research 
that could be studied. First, an important, but classified, thesis topic could involve the 
AJTF PSO algorithm being tested against ISAR returns of actual military targets to verify 
that the motion compensation algorithm can be deployed against such targets of interest. 
Secondly, jet engine modulation (JEM) introduces severe time-variance into an ISAR re- 
turn. Worthwhile research would look at ways of removing the image smearing due to 
JEM or classifying targets based on their JEM signature. A good reference for this re- 
search topic is [15]. 
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APPENDIX MATLAB CODE 


Listed in this appendix is the Matlab code used in this thesis. The order of listing 
is the GA code, followed by the PSO code and the 3D motion detection code. The code 
requires MatLab 6.5 or later and the signal processing toolbox to run. 


As GA CODE 


function [fD1, f{D2, FS, FST] = GA_Image Focus(G, R1, R2,m_orderl,m_order2,num_pop,num_gen,output) 
Written by Capt Wade Brinkman, Canadian Forces 

% Oct 9, 2004 

% Modified: 10 Nov 2004 

% As partial fulfillment of the MSEE Degree at NPS 

% Original DRDC version written by: Jonathan Waisman 

% This version used the exhaustive search technique as described in 

% thesis - see reference page for DRDC paper 


° 
ax 


(=) 


% This algorithm will focus ISAR images using a Genetic Algorithm 
% application of the AJTF algorithm. 


(=) 


% Input: 

% G- signal (2D MXN array) 

% RI -a range bin that will be used for translational motion comp ( 1 to M ) 
% R2-arange bin that will be used for rotational motion comp (1 to M); 
%  m_orderl - order of motion in R1 

% m_order2 - order of motion in R2. 

% num pop - the number in the GA population 

% mnum_ gen - number of generations in the run 

% output - 1 - output TF Transforms 

% 0 - supresses the TF Transforms 

% Output: 

% {D1 - the doppler coefficients for the trans motion comp 

% {D2 - the doppler coefficients for the rotational motion comp 

% FS - final focused signal 

% FST - signal after translational motion comp 


format short; 
[M,N]=size(G); %M -range bins N - pulses 


fD1 = abs(fft(G(R1,:))); 
fD1 = find(fD1 == max(fD1)); 


fD2 = abs(fft(G(R2,:))); 
fD2 = find(fD2 == max(fD2)); 


“find the factorials for 1..m_orderl & two 
for k = 1:m_ order]; 

Orderl_fac(k) = factorial(k); 
end 


for k = 1:m_order2; 
Order2_fac(k) = factorial(k); 
end 


“look at the uncompensated signal; 

if max(size(G)) == 256; 
srdi(G,R1,R2,1,256,1,64,1); colormap(‘hot’); 

else 
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srdi(G,R1,R2,round(N/2)-20,round(N/2)+20,17,31,2); 


end 
*onow we can look at the TF of the two range bins 
if output ==1; 

figure; 


tfr_x1 =tfrsp(G(R1,:)'); “ostft 

tfr_x1 =[tfr_x1(N/2+1:N,:); tfr_x1(1:N/2,:)]; 
imagesc(abs(tfr_x1)); 

title(["Uncompensated - Range Bin # ';num2str(R1)]); 
colormap(‘hot'); 


figure; 
tfr_x2 = tfrsp(G(R2,:)'); “ostft 
tfr_x2 = [tfr_x2(N/2+1:N,:); tfr_x2(1:N/2,:)]; 
imagesc(abs(tfr_x2)); 
title(["Uncompensated - Range Bin # ';num2str(R2)]); 
colormap(‘hot'); 

end; 


WAWL%%V%%%%%%0%%%%%%%%%%%%%%%%%%%%%%%0%%%%%0%%%%%0%%%%%0%%%0%%0% 0% 
%'Step 1: Translational motion compensation using the first range bin' 


[fD1 FST] = GA_AJTF(G,R1,m_orderl,num_pop,num_gen,0,Orderl_fac,fD1); 


WWD %%V%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
W%%%% 
%'Step 2 - Display Translational Motion Compensation Results' 
if max(size(G)) == 256; 
srdi(FST,R1,R2,1,256,1,64,1); colormap(‘hot’); 
else 
srdi(FST,R1,R2,round(N/2)-20,round(N/2)+20,17,31,1); 
end 


“display the 2 range bins after translational rotation 


if output == 
X1=FST(R1,:); Yolst range bin 
tfr_xl =tfrsp(X1'); “ostft 
tfr_x1 =[tfr_x1(N/2+1:N,:);tfr_x1(1:N/2,:)]; 
figure; imagesc(abs(tfr_x1)); 
title({"Range Bin ';num2str(R1),' after Translational Motion Compensation']); 
colormap(‘hot'); 


X2 = FST(R2,:); Yolst range bin 
tfr_x2 =tfrsp(X2'); “ostft 
tfr_x2 =[tfr_x2(N/2+1:N,:);tfr_x2(1:N/2,:)]; 
figure; imagesc(abs(tfr_x2)); 
title({"Range Bin ';num2str(R2),' after Translational Motion Compensation']); 
colormap(‘hot'); 
end; 


%'Step 3 - Rotational motion compensation using the second range bin' 
[fD2 FS] = GA_AJTF(FST,R2,m_order2,num_pop,num_gen,1,Order2_fac,fD2); 
%'Step 4 - Display Rotational Motion Compensation Results' 


if output == 1; 
X1=FS(R1,:); Y%lst range bin 
tfr_xl =tfrsp(X1'); “ostft 
tfr_x1 =[tfr_x1(N/2+1:N,:);tfr_x1(1:N/2,:)]; 
figure; imagesc(abs(tfr_x1)); 
title({"Range Bin ';num2str(R1),' after Rotational Motion Compensation']); 
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colormap(‘hot'); 


X2 = FS(R2,:); Yolst range bin 
tfr_x2 =tfrsp(X2'); “ostft 
tfr_x2 =[tfr_x2(N/2+1:N,:);tfr_x2(1:N/2,:)]; 
figure; imagesc(abs(tfr_x2)); 
title({'Range Bin ';num2str(R2),' after Rotational Motion Compensation']); 
colormap(‘hot'); 
end; 


“finally the compensated image: 
if max(size(G)) == 256; 
srdi(FS,R1,R2,1,256,1,64,1); colormap(‘hot'); 
else 
srdi(FS,R1,R2,round(N/2)-20,round(N/2)+20,17,3 1,2); 
end 


function [fD,NS] = GA_AJTF(SIG,R,m_order,num_pop,num_gen,ToR,Order_fac,fD); 


%Capt Brinkman, Canadian Forces 

%Oct 23, 2004 

%Modified: 20 Nov 04 

“this function is main Genetic Algorithm function that calls all the 
%other GA functions 

“The goal of this real valued GA algorithm is to converge very quickly to 
“the doppler parameters of the basis function 


%SIG: the signal to be analyzed; 

%R: range bin 

%m_order: order of motion 

%num_ pop: this is the population of our genetic pool 

%num_ gen: the number of max # of generations to run 

%ToR: 0 - compensating translational motion 

% 1 - compensating rotational motion 

%Order_fac - precalculated factorials needed for the basis function 
%fD - f1 search parameter 


% Output - fD - returns the doppler coefficients for the motion comp 
% NS - motion compensated signal 


[rbins pulses] = size(SIG); 
index = 1; 


%Step 1: Find the initial Population and rank it 

population = I POP(m_order, num_pop,pulses,fD); 
“%creation our initial population of size num_pop 
%with values from -pulses:pulses 

[gen fit] = FITNESS FFT(population,SIG(R,:),ToR,Order_fac); 

“%ofind the fitness of the first population 

mut_gen = num _ gen; 

while num_gen>0 


gen_score(index) = sum(gen_ fit); 

“ranks the overall fitness of the population 
generation_value(index) = sum(gen_fit)./num_pop; 

%sum of the parents raw score - this should increase over time 
best_value(index) = max(gen_fit); 

%this is the first best value found; 


%Step 2: Based on Suitability, Generate Children 
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children = CHILDS(population,gen_fit,gen_score(index)); 


%Step 3: Mutation 
mut_prob = (0.3); %30% mutation rate. 
children = MUTATION(children,mut_prob,pulses,m_ order); 


%Step 4: Determine Suitability of Children: 
[fit child] = FITNESS FFT(children,SIG(R,:),ToR,Order_fac); 


%Step 5: Reinsertion 
“inserts the best of old population and new population into the new 
% population 
[population,gen fit] = REINSERTION(population,children,gen_fit,fit_child); 


num _gen=num gen -l; 
index = index +1; 
end; 


%pick out the best function; 
best = population(find(gen_ fit == max(gen_fit)),:); 
fD = best(1,:); 


“Now the new signal can be calculated 
if ToR == 0; “translational motion 
hp = basis2(fD,Order_fac,m_order,1,pulses,ToR); 
if pulses<2%oturn off for now 
tfr_hp = abs(stft(hp')); 
tfr_hp = [tfr_hp(pulses/2+1:pulses,:); tfr_hp(1:pulses/2,:)]; 
figure;imagesc(abs(tfr_hp)); 


end; 

NS = SIG.*conj(repmat(hp,rbins, 1 )); 
end; 
if ToR == 1; “rotational motion comp 


clear theta; 
[dummy,theta]=basis2(fD,Order_fac,m_order,1,pulses, ToR); 
spacing = ((max(theta)-min(theta))/pulses)*0.999; 
uniform samples = min(theta)+(1:pulses)*spacing; 
NS = interp1q(theta.';SIG.',uniform_samples.').'; Yoreformat the radar data 
%so that is uniformly samples in theta 
end; 


figure; 
plot(generation_value); 
hold on; plot(best_value,'g.-'); 


function [population] = I POP(m_order,num_pop,N,fD); 


“function I_ POP will create the initial population of possible solutions 
%m_order - order of motion 

’num_ pop - max number in the population 

%N - the values of the pop will be between -N:N; 

%{D - is our fl coefficient as determined by the FFT 


%if m_order == 2; 

% population = linspace(-N,N,num_pop)'; 

%elseif m_order>2; 

% population = [linspace(-N,N,num_pop)' (-N + 2*N*rand(num_pop,m_order-2))]; 
%end; 


population = -N + 2*N*rand(num_pop,m_order-1); 
“guess at f2, f3, f4.....m_order as random vars 
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fD_p(1:num_pop,:) = fD; 
population = [fD_p population]; 
“add f1 to the population 


function [gen_fit] = FITNESS FFT(pop, X, ToR,fac); 
format long 

%Capt Brinkman, 16 Oct 2004 

%Modified 20 Jan 2005 

“%this function generates the suitability of each 
“parent in the current generation 


“pop - matrix of the current population 

% - fl, f2, f3 etc depending on the degree; 
%X - the range bin under examination 

%TOoR - 0 - translational 

% - | - rotational 

“%ofac - precalculated factorials 


“what we are doing here is trying to drive the slope 
“of the scatter in the T - F plane to zero. When this 
% is achieved, translational or rotational motion can be 
“removed. 


[num_pop, degree] = size(pop); 
[a pulses] = size(X); “onumber of pulses in the range bin 


if ToR ==0; %compensating translational motion: 
hp_vec = basis2(pop,fac,degree,num_pop,pulses, ToR); 
“%hp_vec is the basis func for the pop 
X = repmat(X,num_pop,1); 
test_result = X .* conj(hp_vec); 
test = (abs(fft(test_result')).%2)';Yofind the power spectral density of pop 
for k = 1:num_pop; 
test_fft = test(k,:); 


max_test = max(test_fft); 
avg PSD = max _test/sum(test_fft); %Percentage of Power under max point 


gen_fitl(k,l)=avg PSD; %max_ sum; 
“when the percentage of energy in the dominant scatterer of the 
%range bin is a maximum, we have our best focus 

end; 


gen_fit = 100.“(gen_fit1); 
gen_value = gen fitl; 


theta = []; 
end; 


if ToR ==1; Y“ocompensting rotational motion 
[dummy theta] = basis2(pop,fac,length(fac),num_pop,pulses,ToR); 
for k = 1:num_pop; 
spacing = ((max(theta(k,:))-min(theta(k,:)))/pulses)*0.999; 
uniform samples = min(theta(k,:))+(1:pulses)*spacing; 
test_sig = interp1q(theta(k,:)',X',uniform_samples')'; 
“reformat the signal so that it is linear sample w/ 
“theta not time. 
test_fft = abs(fft(test_sig)).“2; “ofind the power spectral density of the scatterer rb 
max_test = max(test_fft); 
avg PSD = max_test/sum(test_fft); %Percentage of Power under max point 
gen_fitl(k,1)=avg PSD; %omax_sum; 
“when the percentage of energy in the dominant scatterer of the 
%range bin is a maximum, we have our best focus 
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end; “end for loop 
end; 


gen_ fit = 100.*(gen_fit1); 
“the exponential provides separation between good values 
“%and bad values. The higher the exponential, the greater the 
%GA searches around first global max found. The smaller, the 
“longer it takes the pop to base towards global max and a 
“global search is executed. 

gen_value = gen fitl; 


function [children] = CHILDS(population,gen_fit,gen_ score); 
%Capt Wade Brinkman 
% Canadian Forces - 20 Oct 04 


“this function will select the best parents based on their assessed fitness 
%and breed them in order to create the children; 


%a roulette wheel selection method is used. This is discussed in detail in 
% Goldberg and the GA toolbox manual 


% population - the current GA pop that we are working with 
%gen_ fit - score between | and zero that assesses their fitness 
“%gen_score - the score for the current gen (sum(gen_fit)); 

[N degree]= size(population); 

%Step 1: Assign Percentage of the roulette wheel 


raw_ score = floor((gen_fit./gen_score).*10000)./10000; 
*% converts gen_fit to a percentage truncated at 4 decimals; 

wheel = cumsum(raw_ score); 
“the values progressively increase and this becomes our random 
“%roulette wheel 


%Step 2: find the parents chosen by the roulette wheel 
select = max(wheel)*rand(1,length(gen_fit)); 

“these are our pointers for the entire selection process 
select = repmat(select,N,1); 

“%convert to matrix 
wheel = repmat(wheel,1,length(gen_fit)); 

*convert to matrix for vector 


choice = wheel>=select; “oplaces a 1 where this is true 
“othe index of the first value in each column corresponds to 
%a parent that has been chosen by the random roulette wheel 


for index = 1:N; 
check = find(choice(:,index)==1); 
Pool(index,:)=population(check(1),:); Y“onecessary for version 6 
“check = find(choice(:,index)==1,1,'first'); -works for version 7 
%Pool(index,:)=population(find(choice(:,index)==1,1,'first'),:); 
“this is our list of good parents that will produce good offspring 
end; 


pooll = Pool(1:N/2,:); 
pool2 = Pool(N/2+1:N,:); Yobreak the pool up into two equal matrix 


“othe first one in the population, after first gen best parent, is 
“guaranteed two spots - one in each pool 
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pooll(1,:) = population(],:); 
pool2(1,:) = population(],:); 


index] = randperm(N/2); 
index2= randperm(N/2); 
Yogenerate a random order; 


alpha =[ones(N/2,1) (0 + 1.*rand(N/2,degree-1))]; Yoalpha is our scaling factor 
“alpha = repmat(alpha,1,degree); “oreformat it into a N/2 - degree 
“%atray to allow for vectorization 
child1 = alpha.*pool1 (index 1,:)+(1-alpha).*pool2(index2,:); 
child2 = (1-alpha).*pool1 (index1,:)+(alpha).*pool2(index2,:); 
“% calculate our children using the method described in 
“thesis and Junfei Li and Hao Ling's paper - see references 


children = [child1;child2]; 
*% combine our children array and return. 


function [population,gen_ fit,gen value] = REINSERTION(population,children,gen_fit.... 
fit_child,gen_value,value_child); 


%Capt Wade Brinkman 
“Modified: 20 Nov 04 
“Input 


“population - old population 

“children - the children of the current population 
%gen_fit - fitness of each element in population 
“%ofit_child - fitness of the children 

%Output 

*%population - new population 

%gen_ fit - new generation fitness 

%gen_value - new generation value 

% This file will create a new population based on fitness 


Y%ofirst step is to combine the two vectors 
temp_ pop = [population;children]; 
temp_fit = [gen_fit;fit_child]; 


for k=1:length(population)/2; 
max_fit = max(temp fit); “ofind the maximum fitness 
max_index = find(temp_ fit == max_fit); %omaximum index 
max_index = max_index(1); “ensures one value is found 


population(k,:) = temp _pop(max_index,:); 
gen_fit(k,:) =temp_fit(max_index); 


temp_fit(max_index(1)) = 0; 
end 


index_new = find(temp_ fit~=0); 
new_fit(1:length(index_new))=temp_fit(index_new); 
new_pop(1:length(index_new),:) = temp_pop(index_new,:); 


for k = (length(population)/2+1 ):length(population); 
index_rand = round(1+(length(population)-1)*rand); 
population(k,:) = new_pop(index_rand,:); 
gen_fit(k,:) = new_fit(index_rand); 

end; 


“%othis 1s our new population with our new measures of fitness; 


function [show_image] = srdi(sig,rb1,rb2,xmin,xmax,ymin, ymax, itype); 
%S1 = srdi(FS,28,19,505,520,17,31,2); will work for fsmall data set 
rbinl = sig(rb1,:); 

rbin2 = sig(rb2,:); 
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show_image = rdi(sig); 


subplot(2,2,[1 3]); 
imagesc(show_image); 

ylabel(‘Down Range - Range Bins’); 
xlabel('Cross Range - Doppler Bins’); 


axis([xmin xmax ymin ymax ]); 
if itype == 0; 

title‘ Uncompensated Image'); 
elseif itype == 1; 

title(Image After Translational Motion Compensation’); 
elseif itype == 2; 

title('Final Image after Rotational Motion Compensation’); 
end 


f rbinl = fliplr(abs(fftshift(fft(rbin 1 ))).%2); 
f rbin2 = fliplr(abs(fftshift(fft(rbin2))).%2); 
n= 2:(length(rbin1)+1); 


subplot(2,2,2); 

plot(n,f_rbin1,'r.-'); 

title(["Power Spectrum Range Bin ';num2str(rb1)]); 
axis([xmin xmax min(f_rbin1) max(f_rbin1)]); 
ylabel(‘Amplitude'); 

xlabel('Cross Range - Doppler Bins’); 


subplot(2,2,4); 

plot(n,f_rbin2,'r.-'); 

title({"Power SpectrumRange Bin ',;num2str(rb2)]); 
axis([xmin xmax min(f_rbin2) max(f_rbin2)]); 
ylabel(‘Amplitude’); 

xlabel('Cross Range - Doppler Bins'); 


function [hp,theta] = basis2(pop,fac,degree,xnum_pop,pulses, ToR) 

%Capt Wade Brinkman 

%16 Oct 04 

“Modified 20 Jan 05 

“this function will calculate the basis function for our current population 
“%the calculation will be done in parallel to speed up comp time 

“this version removes the for loops which cuts down the run time by half 
“this 1s important since we are trying to drive the run-time of the GA and Swarm to 
“almost zero 

“Input: 

%pop - current population 

“%ofac - the factorials used in the calculation 

“degree - the degree of motion 

’num_pop - the number in current population 

“pulses - # of pulse in cross range 

%ToR - Translational or rotational 0 - theta = []; 1 - hp = []; 

%Output: 

“hp - the basis function of each mbr of pop 

“theta - the phases associated with each mbr of pop 

hp = zeros(num_ pop, pulses); 


t = ((1:pulses)./pulses)'; Yot = [(1 2 3 4... pulses)./pulses]' - column vector 
t_temp = t; 
for 1 = 2:degree 
t = [t (t_temp).”1]; 
%t’*1 forms a column of t, t2 forms the next...used to vectorize 
end 


pop_fac = pop./(repmat(fac.num_pop,1)); “opredividing the population by its factorials 
theta = (pop_fac(1:num_pop,:)*t(1:pulses,:)'); 
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if ToR == 1; 
hp = []; 
else 
hp(1:num_pop,1:pulses) = exp(theta.*(-j*2*p1)); 
theta = []; 
end; 
“finally, this calculates the basis function for the entire population. 
function f= rdi(s) 
% generates a range doppler image of the signal s 


f = fft(s'); 
[n,m] = size(f); 
f= [f(round(n/2)+1:n,:); {1 :round(n/2),:)|; % equivalent to fftshift 


f = abs(f/1); 
f =( f./max(max(f))); 


figure; imagesc(abs(f)); 
colormap(‘jet'); 


B. PSO CODE 


function [FS fDT fDR] = Swarm Image Focus 3D(G, R1,D_T, RFocus,D_R,num_particles,cycles) 


% Written by Capt Wade Brinkman, Canadian Forces 

% Nov 13, 2004 

% Modified: 12 Mar 2005 

% As partial fulfillment of the MSEE Degree at NPS 

% Original DRDC version written by: Jonathan Waisman 

% The DRDC version used the exhaustive search technique as described in 
% thesis - see reference page for DRDC paper 

% This algorithm will focuse ISAR images using a Swarm Algorithm 

% application of the AJTF algorithm. 

% This algorthm uses a Swarm Algorithm as described in Jacob Robinson 
% and Daniel Boeringer papers - see reference page - with modifications 
% to fit the AJTF algorithm 

% This algorithm uses several scatterers from the same range bin for 

% rotational motion compensation. 

% Input: 

% G- signal (2D MXN array) 

% RI -arange bin that will be used for translational motion comp 

% D_T - degree of translational motion compensation 

% RFocus - [Range Bin, # Scatterers, Range Bin, # Scatterers,....] 

% list of range bin and number of scatterers for rotational motion 

% compensation 

% D_R - degree of rotation motion compensation 

% num_ particles - the number of particles in the Swarm - identical to 

% population 

% cycles - number of times the Swarm's position are updated - identical 
% to number of generations in GA 


% Output: 

% {DT - the doppler coefficients for the translational motion 

% compensation 

% {DR - the doppler coefficients for the rotational motion compensation 
% FS - final focused signal 


format short; 

[M,N]=size(G); %M -range bins N - pulses 

N_Rot = length(RFocus); 

if mod(N_Rot,2)~=0; 
disp(‘Rotational Range Bin Matrix incorrectly formatted’); 
FS = G;fDT = []; {DR = []; 
return 

end 
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N_RSB=N _ Rot/2; “number of range bins for rot mot comp 
RFocus = reshape(RFocus,2,N_RSB); 


for k = 1:D_T;%find the factorials 
Orderl_fac(k) = factorial(k); 


end 
fork = 1:D_R; 

Order2_fac(k) = factorial(k); 
end 


WAL %V%%%%%%0%%%%%0%%%%%0%%%%%%%%%%%%%0%%%%%0%%%%%0%%%%%0%%%%%% 0% 
%'Step 1: Translational motion compensation using the first range bin' 

%For this step, the Swarm _AJTF_ 3D will return the Doppler Frequencies, but not 

“focus the image 


[fDT] = Swarm _ AJTF 3D(G,R1,D_T,num_particles,cycles,0,Order1_ fac); 
FST = tr_focus_3D(G,fDT,0); 
WWMM %%%%%0%%%%VM%%%%%%%%%%%%%%%%MMM%%%%%%%%%%%%% 


FST_t=FST; 
%'Step 2 - Rotational motion compensation using the Range Bins 
’%and number of scatterers found in RFocus 
r fDR=1; “index for the rotational doppler 
fork = 1:N_RSB 
for g = 1:RFocus(2,k); 
[fDR(r_fDR,:)]=Swarm_AJTF 3D(FST_t,RFocus(1,k),D R,num_particles,cycles,1,Order2_ fac); 
if RFocus(2,k)>1&&g<RFocus(2,k) 
X_t=fft(FST_t(RFocus(1,k),:)); 
index = f{DR(r_fDR,1)+(-3:3); 
for kk = 1:length(index); 
if index(kk)<1 X_t(N+index(kk)) =0; 
elseif index(kk)>N X_t(index(kk)-N)=0; “oremove used scatterer from Range Bin 
else X_t(index(kk)) =0; end; 


end; 
FST_t(RFocus(1,k),:) = ifft(X_t); 
end; 
r fDR=r fDR+1; 
end 


end; 


FD = mean(fDR, 1); 
FS = tr_focus_3D(FST,FD,1); 
function [fD] = Swarm _AJTF_ 3D(SIG,R,m_order,num_particles,cycles,ToR,Order_fac,fD); 


%Capt Brinkman, Canadian Forces 

%Nov 13, 2004 

“Modified 20 Jan 05 

“%this function is main Swarm Algorithm function that calls all the 
“other Swarm functions. It is called once for translational and once for 
“%rotational motion compensation. 

%The goal of this real valued Swarm algorithm is to coverge very quickly to 
%the doppler parameters of the basis function 

%this is a line search - it will find fl, fix f1, find f2, fix f2 find f3 

“etc. Necessary since solution space grows by N*m_ order 

%SIG: the signal to be analysed; 

%R: range bin 

%m_order: order of motion 

*%num_ particles: this is the number of partilces in our Swarm 

%ocycles: max # of generations to run 

%ToR: 0 - compensating translational motion 


% 1 - compensating rotational motion 
%Order_fac - precalculated factorials needed for the basis function 
%fD = fl 


% Output - fD - returns the doppler coefficients for the motion comp 
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fD = abs(fft(SIG(R,:))); f{D = find(fD == max(fD)); 
[rbins pulses] = size(SIG); 

num_cycles = cycles; 

GFit_ Old = 0; 


current f= 2; Yocurrent parameter being searched 
Swarm=zeros(num_particles,1); 
Swarm(1:num_particles,1)=[fD]; 


while current_f<=m_order 


%Step 1: Find the initial Swarm and evaluate fitness 
[Swarm,velocity] = 1 Swarm(Swarm,m_ order, num_particles,pulses,current_f); 
“creation our initial swarm and velocity 
% with values from -pulses:pulses and intial velocities from 
%-pulses:pulses 
“%eval the fitness 
[part fit] = FITNESS FFT(Swarm,SIG(R,:),ToR,Order_fac(1:current_f),current_f); 
LBest = Swarm; “these are the local 'bests' found by each particle. For the first cycle, by default it is 
%equivalent to the first Swarm 
LFit = part fit; %LFit is the fitness of the local best 
GFit = max(LFit); Yothis is the maximum global fitness that was found 
GBest = LBest(find(LFit==max(LFit)),:); “othis is current global best...... 
GFit = GFit(1); GBest = GBest(1,:); 
index = 1; 
part_score(index) = sum(LFit)./num_ particles; 
*%ranks the overall fitness of the population 
best_ value(index) = GFit; 
“this is the current global maximum; 


best_count = cycles; 
o Best = GFit; “these two variables are used to cause the algorithm to exit if no improvement has been made after x - cycles 
“figure; 
while cycles>0 
“hold on; plot(index,LBest(:,2),'g+');hold on;plot(Gindex,GBest(1,2),"b>');hold on; plot(index,Swarm(:,2),'c.'); 
cycles = cycles -1; 
index = index +1; 


%Step 2: Update velocities: 
[velocity]=v_update(velocity,_LBest,GBest,swarm,num_particles,current_f,pulses); 


%Step 3: Update positions 
[Swarm,velocity]=p_update(Swarm,velocity,num_particles,pulses); 


%Step 4: determine fitness of the Swarm 
[part_fit] = FITNESS FFT(Swarm,SIG(R,:),ToR,Order_fac(1:current_f),current_f); 


%Step 5: determine LBest and GBest 
[LBest, LFit, GBest, GFit] = swap_best(LBest, LFit, GBest, GFit, Swarm,part_fit,current_f); 


part_score(index) = sum(LFit)./num_ particles; 
“ranks the overall fitness of the population 
best_value(index) = GFit; 
“this is the current global maximum; 


end; 


if GFit_Old >= GFit 
GBest(current_f)=0; “if the Fitness of this degree is less than the previous, set this parameter 
“equal to zero. 
else 
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GFit_Old = GFit; 
end 


current_f=current_f+1; 
cycles = num_cycles; 
Swarm=repmat(GBest,num_particles, 1); 


end; 


fD = GBest; 


function [Swarm, velocity] = 1 Swarm(Swarm,m_order,num_particles,N,current_f); 
%Capt Wade Brinkman, Canadian Forces 

%14 Nov 04 

%Modified 12 Jan 05 

“function I_ Swarm will ceate the initial particles of possible solutions 

“%and assign velocities random in both mag and direction to each particle 


%Swarm - current Swarm - will be expanded 

%m_order - order of motion 

’num_particles- max number in the particles in the Swarm 
%N - the values of the pop will be between -N:N; 

%{D - is our fl coefficient as determined by the FFT 


Swarm =[Swarm (-N + 2*N*rand(num_particles,1))|; 
“guess at f2, f3, f4.....m_order as random vars 

velocity = (2*N)*rand(num_particles, 1); 
%magnitude of the velocites 

v_dir = rand(num_particles,1)-0.5; 
“%direction of the particle +ve or -ve 

index! = find(v_dir>=0); 

index2 = find(v_dir<0); 

v_dir(index1)=1; 

v_dir(index2)=-1; 


velocity = velocity.*v_dir; 
velocity = [zeros(num_particles,current_f-1) velocity]; 
“the f1 velocity will always be fixed at zero. 


function [velocity] = v_update(velocity,LBest,GBest,Swarm,num_particles,m_order,pulses); 
%written by Capt Wade Brinkman, Canadian Forces 

%Created 14 Nov 2004 

% This function will update the velocities based on the particles 
*%current velocity and the pull of the the Local Best and the GBest 
“Input: 

“velocity - current velocity vectors 

%LBest - matrix of local best discoveries 

%GBest - the best discovery 

% Swarm - current search location of the particles 

*%num_particles - number of particles in the swarm 

%m_ order - order of motion 

“pulses - number of pulses - also corresponds to max velocity 
%Output: 

%velocity - updated velocity vector 


K(1:num_particles,1:m_order)=0.729; “constriction factor 
rhol = 2.05; 

rho2 = 2.05; 

r_vecl =rand(num_particles,m_order); 

r_vec2 = rand(num_particles,m_order); 

Global Best=repmat(GBest,num_particles,1); 
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velocity = K.*(velocity + 2.05.*r_vecl.*(LBest-Swarm)+2.05.*r_vec2.*(Global_Best-Swarm)); 
“this sets the velocity vector for the particles 
“Now, check to make sure that no particle exceeded the max velocity; 
[R C V] = find(velocity>2*pulses); 
if length(V)>0; 

for k=1:length(R) 

velocity(R(k),C(k)) = 2*pulses; 

end 

end 


[R C V] = find(velocity<-2*pulses); 
if length(R)>0 
for k=1:length(R) 
velocity(R(k),C(k)) = -2*pulses; 
end 
end 


function [Swarm,velocity]=p_update(Swarm,velocity,num_particles,pulses); 
%Programmed by Capt Wade Brinkman, Canadian Forces 

%14 Nov 2004 

“based on the new velocities, the positions will be updated. Any particle 
“outside the solution space of +/- pulses will have its velocity in that 

% dimension zeroed 

“Input: 

%Swarm: The particle Swarm 

%velocity: Particles Current velocity 

*%num_ particles: Number of particles in the Swarm 

“pulses: number of pulses in signal which also defines the bounds 

% of the solution space 


Swarm = Swarm + velocity; 
[R C V] = find(Swarm>pulses); 


if length(V)>0; “ot+ve wall 
for k = 1:length(V); 
velocity(R(k),C(k))=-pulses*0.1*rand; Yorebounding wall 
Swarm(R(k),C(k))=pulses; 
end 
end 


[R C V] = find(Swarm<-pulses); 


if length(V)>0; %-ve wall 
for k = 1:length(V); 
velocity(R(k),C(k))=pulses*0.1*rand; “orebounding wall 
Swarm(R(k),C(k))=-pulses; 
end 
end; 


function [LBest, LFit, GBest, GFit] = swap_best(LBest, LFit, GBest, GFit, Swarm,part_fit,m_order); 
% Written by Captain Wade Brinkman, Canadian Forces 

%14 Nov 2004 

% This function will update the local and global bests as required 
%oInputs: 

%LBest - matrix of local best parameters 

%LFit - a column vector of the fitness of the local bests 

%GBest - is the parameters of the global best 

%GFit - fitness of the global best; 

% Swarm - current parameters of the particle swarm 

%part fit - current Swarm's fitness 

%m_order - order of motion 

%Output: LBest, LFit, GBest, GFit - updated parameters 


Max_Swarm_fit = max(part_fit); othe best of the current particles 
Max Swarm = Swarm(find(part_ fit ==Max_ Swarm _fit),:); 
if Max_Swarm_fit>GFit; 
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GFit = Max Swarm _fit(1); 
GBest = Max_Swarm(1,:); 
end; 


vec] = LFit>=part_fit; “ois current LBest better than current location 
vec2 = part_fit>LFit; “is current location better than LBest 


LFit = LFit.*vecl+part_fit.*vec2; Yoreplaces the better current positions and LFit 
LBest = LBest.*repmat(vec1,1,m_order)+Swarm.*repmat(vec2,1,m_order); “oplaces the better position in Local Best 


function [Focused_ Signal] = tr_ focus 3D(signal,fD,ToR); 
%Capt Wade Brinkman 

%15 Nov 2004 

“this is a utility function tha will perform translational 
“%or rotational motion compensation on the signal 


[rbins pulses] = size(signal); 
for k = 1:length(fD); 

o_fac(k) = factorial(k); 
end 


[hp theta] = basis2(fD,o_fac,length(fD),1,pulses, oR); 


if ToR = 0 
Focused_ Signal = signal.*conj(repmat(hp,rbins, | )); 
else 


spacing = (max(theta)-min(theta))/pulses*0.99999; 

uniform samples = min(theta)+(1:pulses)*spacing; 

Focused_ Signal = interp1q(theta.',signal.';uniform_samples.').'; 
end; 


c. 3D MOTION DETECTION CODE 


function [NL]=MDetect_3D(namedata,num_pulses,load_data,save_data,output,Num_S,ord1,ord2,cmovie); 
“function [NL]=MDetect_3D(‘data',2*11,1,0,1,2,2,2,0); 

%Capt Wade Brinkman 

%Created 6 Feb 2005 

“Modified 12 March 2005 

% This is a program that will detect the presense of 3D motion 

“%It uses the EAJTF Algorithm {as described by Thayaparan's unpublished work} 
*%with the Swarm algorithm as its search engine 

% Linearity of Phases is the Average Linearity of Phases Method 

% {as described by Thayaparan's unpublishedd work} 

%INPUT 

% namedata - the name of the data set {without the imaging interval 

% number} 

% num_pulses - number of ISAR pulses in the cross-range 

% - either 2“10 or 2“11 

% load_data - if 1 - fD and theta already exisit 

% save data - save the fD and theta 

% output - output all images 

% Num_S - number of Scatterers in first range bin {default 1 1n second} 
% ordl - order of motion comp - Translational 

% ord2 - order of motion comp - Rotational 

% cmovie - create a movie 

% OUTPUT 

% NL - Degree of Non-Linearity 


fflag = 1; 
num_images = floor(60000/num_pulses); 


“first set the directory of choices 
if (num_pulses == 2“11)&&(ord1==2); 
cd LP2048 2D; 
directory = cd; cd ..; “sets our active directory where fD and theta are stored 
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rb_file='rb focus 3D'; “this is the file that contains the range bins used 

temp_load = cat(2,directory,'/',rb_ file); 

load(temp_ load); “%the rb_ focus 3D contains the range bins for 3D motion detection 
rb focus =rb_ focus 3D; clear rb focus 3D; 


elseif (num_ pulses == 2%11)&&(ord1==4); 


cd LP2048 4D; 

directory = cd; cd ..; %sets our active directory where fD and theta are stored 
rb_file='rb focus 3D'; “this is the file that contains the range bins used 

temp_load = cat(2,directory,'/',rb_ file); 

load(temp_ load); “the rb_ focus 3D contains the range bins for 3D motion detection 
rb_ focus =rb_ focus 3D; clear rb focus 3D; 


elseif (num_pulses == 2“10)&&(ord1==2); 


cd LP1024 2D; 

directory = cd; cd ..; %sets our active directory where fD and theta are stored 

rb_file ='rb focus2'; “othis is the file that contains the range bins used 

temp_load = cat(2,directory,'/',rb_ file); 

load(temp_load); “%the rb_ focus 3D contains the range bins for 3D motion detection 
rb_ focus =rb_focus2; clear rb focus?2; 


elseif (num_ pulses == 2%10)&&(ord1==4); 


cd LP1024 4D; 

directory = cd; cd ..; %sets our active directory where fD and theta are stored 
rb_file='rb focus2'; “othis is the file that contains the range bins used 

temp_load = cat(2,directory,'/',rb_ file); 

load(temp_load); “the rb_ focus 3D contains the range bins for 3D motion detection 
rb_focus =rb_focus2; clear rb focus?2; 


else 


disp('Must choose 2“11 or 2“10 pulses or Order of Motion Error equal to 2 or 4'); 
return 


end; 


for k =1:num_images; 
if load_data == 0; “oneed to generate new fDs and phases 
if fflag == 1; fflag = 0; fac = zeros(1,ord2); for g = l:ord2; fac(g) = factorial(g); end; end; “oprecalc factorials 
image_int = cat(2,directory,'/',namedata,int2str(k-1)); Yocreate string containing name of imaging interval 
loadGmage int); “load imaging interval into work space 
[FS fDT fDR] = Swarm Image Focus 3D(dat,rb_focus(k,1),ord1,[rb_ focus(k,2),Num_S,rb_focus(k,3),1 ],... 
ord2,20,100); 
if output==1 
rdi(FS); 
title({'Focused Image Imaging Interval ';num2str(k-1),'- 'jnum2str(num_pulses),' Pulses in Cross Range']); 
xlabel('Cross Range - Doppler Bins’); 
ylabel('Down Range - Range Bins’); 
end 
[dum ord2] = size(fDR); f = [fDR]; [m n] = size(f); 
[hp,theta] = basis2(f,fac,ord2,m,num_pulses, 1); 
elseif load_data == 1; “data already exisits 
doppler = cat(2,directory,'/','fD',int2str(k-1)); 
load(doppler); 
[r c] = size(fD); 
fDT = fDC,:); 
fDR = fD(2:r,:); 
if output == 1; 
image_int = cat(2,directory,'/',namedata,int2str(k-1)); 
load(image_int); 
FST = tr_focus_ 3D(dat,fDT,0); 
FS = tr_focus(FST,mean(fDR, 1),1); 
title({'Focused Image Imaging Interval ';num2str(k-1),'- 'j;num2str(num_pulses),' Pulses in Cross Range']); 
xlabel('Cross Range - Doppler Bins'); 
ylabel('Down Range - Range Bins’); 
end; 
[dum ord2 |=size(f{DR); Y%odetect number of scatterers and order 
phases = cat(2,directory,'/','theta',int2str(k-1)); 
load(phases); 
if fflag == 1; fflag = 0; fac = zeros(1,ord2); for g = 1:ord2; fac(g) = factorial(g); end; end; “oprecalc factorial 
f = [fDR]; [m n] = size(f); 
end; 
[NL(k)] = E_DEG OF NONLINEARITY 2( theta ,f;min(size(theta)),fac,0 ); 
if save_data == 1; 
doppler = cat(2,directory,'/','fD',int2str(k-1)); 
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fD = [fDT;fDR]; 
save(doppler,'fD'); 
phases = cat(2,directory,'/','theta',int2str(k-1)); 
save(phases, 'theta’'); 
end; 
end; 


NL _norm = NL; 
figure; plot(O:(num_images-1),NL_norm,'r+-','linewidth',2); grid on; axis tight 
cutoff line(1:(num_images)) = mean(NL); 
“hold on; plot(0:(num_images-1),cutoff_line,'b-.',"linewidth',2); hold off; 
xlabel(‘Imaging Interval’); 
ylabel(‘Degree of Non-Linearity'); 
title({'Degree of Non-Linearity in Each Imaging Interval"... 
['6 - Point Delta Wing Data Set with ';num2str(num_pulses).... 
"Pulses and Degree of Motion Compensation - ';num2str(ord1)]}); 


%now we can display the good images 
%and make a movie if the option is selected 
[rc v] = find(NL_norm<=.010*mean(NL)); 
for k =1:length(c); 
image _ int = cat(2,directory,'/',namedata,int2str(c(k)-1)); 
load(image_ int); 
doppler = cat(2,directory,'/','fD',int2str(c(k)-1)); 
load(doppler); 
[row col] = size(fD); 
fDT = fD(1,:); 
fDR = fD(2:row,:); 
fD = mean(fDR, 1); 
FST =tr_ focus 3D(dat,fDT,0); 
FS = tr_focus(FST,fD, 1); 
image = rdi2(FS); 
bound = image>.3; 
[r2 c2 v2] = find(bound==1); 
center_r = round((max(r2)+min(12))/2); 
center_c = round((max(c2)+min(c2))/2); 
axis([center_c-15 center_ct+15 8 38]); 
if cmovie ~=1; title(['Imaging Interval: 'snum2str(c(k)-1)]); end; 
if cmovie ==1; 
axis off; 
pause(2); 
g 04 NL(k) = getframe(gcf); 
end 
end 


if cmovie == 1; save g 04 NL g 04 NL; end; 


%now we are going to plot out the three best and three worst imaging 
Yointervals 
NL t=NL; 
for k = 1:3 “find the three best 
temp = find(NL_t == min(NL t)); 
best_image(k) = temp-1; NL_t(temp) = NaN; 
end 


for k = 1:3 %find the three worst 
temp = find(NL_t == max(NL t)); 
worst_image(k) = temp-1; NL_t(temp) = NaN; 
end 


“%now plot out the 3 best and 3 worst along with their non-linearity plots 


fork = 1:3 
figure; 
image _ int = cat(2,directory,'/',snamedata,int2str(best_image(k))); 
load(image_int); 
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phases = cat(2,directory,'/',‘theta',int2str(best_image(k))); 
load(phases); 
doppler = cat(2,directory,'/’,'fD',int2str(best_image(k))); 
load(doppler); 
[row col] = size(fD); 
fDT = fD(1,:); 
fDR = fD(2:row,:); f(D = mean(fDR, 1); 
FST = tr_focus_ 3D(dat,fDT,0); 
FS = tr_focus_ 3D(FST,fD,1); 
im_int = rdi2(FS); 
subplot(2,2, 1); 
imagesc(im_ int); 
title({'Imaging Interval ',;num2str(best_image(k)),' Low Non - Linearity']); 
bound = im_int>.3; 
[r2 c2 v2] = find(bound==1); 
center_r = round((max(r2)+min(r2))/2); 
center_c = round((max(c2)+min(c2))/2); 
axis([center_c-15 center_ct+15 8 38]); 
[test] = E DEG OF NONLINEARITY 2( theta ,fDR,min(size(theta)),fac,1 ); 
end; 


fork = 1:3 
figure; 
image int = cat(2,directory,'/',snamedata,int2str(worst_image(k))); 
load(image_int); 
phases = cat(2,directory,'/,‘theta',int2str(worst_image(k))); 
load(phases); 
doppler = cat(2,directory,'/','fD',int2str(worst_image(k))); 
load(doppler); 
[row col] = size(fD); 
fDT = fD(1,:); 
fDR = fD(2:row,:); f(D = mean(fDR, 1); 
FST = tr_focus 3D(dat,fDT,0); 
FS = tr_focus_ 3D(FST,fD,1); 
im_int = rdi2(FS); 
subplot(2,2, 1); 
imagesc(im_ int); 
title(['Imaging Interval ';num2str(worst_image(k)),' High Non - Linearity']); 
bound = im_int>.3; 
[r2 c2 v2] = find(bound==1); 
center_r = round((max(r2)+min(r2))/2); 
center_c = round((max(c2)+min(c2))/2); 
axis([center_c-15 center_ct+15 8 38]); 
[test] = E DEG OF NONLINEARITY 2( theta ,fDR,min(size(theta)),fac,1 ); 
end; 
function [AVGN] =E DEG OF NONLINEARITY2(P, fD , L ,fac,OP) 
% AVGN =E DEG OF NONLINEARITY2(P, fD,L) - 
“% Finds the presence of 3D motion in each imaging interval. And returns 
% the following. 
% 
% INPUT(S): 
% P - the phase function of an imaging interval 
% {D - doppler coefficients for this phase function (P can be generated 
% from fD) 
% L - number of point scatterers chosen for analysis 
% 
% OUTPUT(S): 
% AVGN - degree of 3D motion for the imaging interval that defined the 
% phase function P. Note that AVGN is the average degree of 


% non-linearity (or in other words the degree of 3D motion) for the non-linearities 
% for every pair of phase functions defined in P. 
% 


% Author : Jonathan Waisman 

% Modified by Capt Wade Brinkman, Canadian Forces 
% First Modified: 9 Feb 05 

% Last Modified: 29 Mar 05 


[H , N] =size( P ); 
[B , ORDER] = size( fD ); 
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Ny = zeros( 1, L ); 
AVGN = 0; 
AProj = 0; 


“%ideal doppler coef will be average 

i f{D = mean(fD, 1); 

“%ideal phase function 

[hp 1. theta]=basis2(1_fD,fac,ORDER,1,N,1); 

aa = max(max([1i_theta;P])); 

fori=1:L 
% generate the non-linearity between the IdealP and P(,:) 
[Ny(i )] =E_DEG OF NONLINEARITY 1(1_ theta , PG,:),OP,i+1 ); 
% sum the non-linearities for averaging 
AVGN = AVGN + Ni(1 ); 

end 


AVGN = AVGN /L; 

function | N12 Proj] = E DEG OF NONLINEARITY1(i P, r_P,OP,index ) 
% | N12 |= DEG OF NONLINEARITY1(i P,r P ) - calculates the degree of 
% non-linearity between 1 P andr P using a least squares approach. 

% 

% INPUT(S): 

% r_P,i P - 2 phase functions 

% 

% OUTPUT(S): 

% N12 - degree of nonlinearity between these 2 particular phase functions 

% 

% Author : Jonathan Waisman 

% Modified by Capt Wade Brinkman 

% Last Modified: 29 March 05 

% 1 P(t)=a *r P(t) +n(t) is the relation of 2 point scatterers. 


p=polyfita Pr P,1); 
p = flipIr( p ); 


N = length(i_P ); 
k = length( p ); 
func = zeros( 1 , N); 
t = linspace(1,max(i_ P),N); 
for x = 1:N 
for1= 1:k 
func( x ) = func( x )+ p(i) * (i _P(x)*%(1i-1)); 
end 
end 


% p =[p(1),p(2)]| where least squares approximation func(t) = p(1) + p(2)*t 
% but func is the best fit linear approximation to 1_P, and thus: 

% the degree of nonlinearity is the difference in the plot of 

% plot(i_P,func) and plot(i_P,r_P); 


N12 =0; 


P=r_P: 
FUNC = func; 
for x = 1:N 
N12 =N12 + abs(P(x)-FUNC(x))./abs(FUNC(x)); 
end 


if OP == 
subplot(2,2,index); hold on; 
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plotG_P,func,'g.-','linewidth',2);hold on; 
plot(@_P,r_P,'r-.','linewidth',2); 
xlabel(‘\Theta_i1 d_e_a_1'); 
ylabel({['\Theta_';num2str(index-1),'(t)'];' as a Function of \Theta i de a I'}); 
grid on; axis tight 
end; 


“to see the degree of non-linearity, enter this line into the command line: 
“figure; plot(i_P,func,'g.-','linewidth',2);hold on; 
“%plot(i_P,r_P,'r.-','linewidth’,2); 


% The degree of 3D motion as the deviation from a linear 

% relationship between theta and phi over the dwell intervall 

% is directly related to N by aconstant Beta. Thus the degree of 
% nonlinearity will give a degree of 3D motion! 
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